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Polyimide Composites have received considerable attention due to their present and potential 
future application, owing to their small structural weight and higher glass transition 
temperature, Tg. In this work the formation of cracks from smooth surfaces of brittle materials 
was studied. 

It is known from the previously published research that polyimide composite plates are 
sensitive to thermal aging that manifest itself by the degradation of the effective mechanical 
properties and weight loss due to oxidation processes. Microscopical examination revealed 
generation of the microcracks perpendicular to fiber direction and voids in polymer matrix 
material. 

The microcracking generation mechanisms depend on the environmental conditions such as 
time and temperature but also on the load history for a given material. Microstructure of the 
material itself, fiber volume ratio, type of reinforcement, fiber distribution and fiber alignment 
definitely influence the aging and microcracking processes. 

As the first approximation of the aging process the crack formation in isotropic material can be 
studied due to non uniform strain distribution caused by aging phenomenon. Cracking is the 
main failure mode for brittle materials when the stress is tension dominant. The starting point 
of cracking can be determined by various strength criteria. However, the strength criteria does 
not contain any information that dictates the behavior of the following cracking in the material. 

In the following section a copy the Ph.D dissertation of the student Anping Hong supported by 
this NASA project is included as a report of the research efforts. This work can be expanded 
in the future to address the cracking in the orthotropic materials. The theoretical results should 
be first examined using aging of the resin plate that can be modeled by the work done in the 
presented studies. „ 
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ABSTRACT 


Sudden appearance of cracks on a macroscopically smooth -surface of brittle 
materials due to cooling or drying shrinkage is a phenomenon related to many 
engineering problems. Although conventional strength theories can be used to predict 
the necessary condition for crack appearance, they are unable to predict crack spacing 
and depth. On the other hand, fracture mechanics theory can only study the behavior of 
existing cracks. The theory of crack initiation can be summarized into three conditions, 
which is a combination of a strength criterion and laws of energy conservation, the 
average crack spacing and depth can thus be determined. The problem of crack initiation 
from the surface of an elastic half plane is solved and compares quite well with available 
experimental evidence. 

The theory of crack initiation is also applied to concrete pavements. The 
influence of cracking is modeled by the additional compliance according to Okamura’s 
method. The theoretical prediction by this structural mechanics type of model correlates 
very well with the field observation. The model may serve as a theoretical foundation for 
future pavement joint design. 

The initiation of interactive cracks of quasi-brittle material is studied based on a 
theory of cohesive crack model. These cracks may grow simultaneously, or some of 
them may close during certain stages. The concept of crack unloading of cohesive crack 
model is proposed. The critical behavior ( crack bifurcation, maximum loads) of the 
cohesive crack model are characterized by rate equations. The post-critical behavior of 
crack initiation is also studied. 
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CHAPTER I 


INTRODUCTION 

Formation of initial macroscopic cracks from smooth surfaces of brittle 
materials can be frequently observed. For instance, in large blocks of concrete, nearly 
periodic cracks can appear suddenly due to drying shrinkage or hydration heat. In the 
case of reinforced concrete beams, cracks of a certain spacing form at the tensile face 
when the applied load becomes large enough. In drying lake beds and mud flats, cracks 
with a honeycomb pattern emerge after some period of drying. When a floating sea ice 
plate is subjected to a vertical load, star-shaped cracks of finite length radiate suddenly 
from the loaded area when the applied load reaches a certain level. Highway pavement 
develops cracks each year due to various environmental changes such as temperature 
and moisture. In all these situations, cracks of macroscopic sizes form suddenly from a 
smooth surface. 

Cracking is the main failure mode for brittle materials when the stress state is 
tension dominant. The starting point of cracking can be marked by various strength 
criteria depending on the specific material under specific environment. In the simplest 
case, one can use the condition that the maximum tensile stress should not exceed the 
tensile strength of the material. However, the strength criteria does not contain any 
information that dictates the behavior of following cracking in the material. Engineers 
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used to avoid the cracking by reducing the maximum stress in their designs. However, 
cracking in some cases is either inevitable or uneconomical. For instance, the asphalt 
concrete pavement in northern America and Canada develops crack each and every 
winter, it cannot be avoided. It would not be cost effective to design a reinforced 
concrete beam prohibiting cracking on its tension side. Sometimes, it may become 
crucial to know what will happen after cracking occurs. For instance, due to safety 
concern, it is very important to know the behavior of nuclear reactor structure after 
cracking has occurred. Obviously, the conventional strength criteria are not sufficient 
to answer these questions. 

Cracking has been studied by the theory of fracture mechanics. During the 
recent half century since the practical applications pioneered by Irwin, fracture 
mechanics has been significantly developed and enriched in terms of practical 
applicability and theoretical understanding. Although the original intention of the 
fracture mechanics, initiated by Griffith (1924), was to understand the strength of 
brittle materials and its relation with the inevitable defects in the material, the later 
development was more focused on the problems of crack growth. As a result, a new 
type of material strength, called the toughness of the material (or fracture energy), was 
brought into the light and became the subject of intensive research. Fracture toughness 
of the material determines whether a crack of a given size will propagate under a given 
load. Fracture mechanics works where the conventional strength criteria fail, because 
at crack tip the stress is infinite so the conventional strength is inapplicable. However, 
fracture mechanics does not encompass the territory of the conventional strength 
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criteria, either. If there are no pre-existing cracks in material (or the pre-existing cracks 
are very small in sizes), fracture mechanics would predict unreasonably high failure 
load. Mathematically this is because the energy release rate for a very short crack 
approaches zero like a linear function of crack length (so the stress intensity factor due 
to a unit load approaches zero like the square root of the crack length). 

It is, therefore, seen that the theory of fracture mechanics and the strength 
criteria have their own domains of applicability that do not overlap. The phenomenon 
of crack initiation lies partially in both of these two domains, but not entirely in either 
one. This is because the strength theory alone cannot describe the cracking behavior, 
but fracture mechanics cannot be adequately employed either where there is no pre- 
existing cracks. A new theory that combines these two different aspects of material 
strengths must be found to study the crack initiation theory. This kind of theory does 
not seem to have been systematically studied before. 

From the practical point of view, the pattern of initial cracking can be important 
for many engineering problems. Sometimes sparse cracking patterns are preferred, 
sometimes dense but small cracking patterns are more desirable. For the former case 
one can mention pavement cracking. The longer crack spacing is preferred because 
one can use less joints to control cracking at less cost. For the later case, one can 
mention cracking in the reinforced concrete. To prevent corrosion of reinforcement, 
the cracks must remain hairline thin, and so we need the crack spacing to remain 
sufficiently small. 
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In all these cases, it is important to understand the mechanism of crack initiation 
so that we can know what can be done to control cracking. The main focus of this 
dissertation is to establish the formulation of the crack initiation theory, solve the crack 
initiation problems of different geometric structures as well as different material 
properties. It is hoped that these new knowledge can shed new light on many 
engineering problems in which cracking is inevitable. 

Some basic idea of crack initiation theory for the case of vertical penetration of 
a floating sea ice plate was contained in the study by Li and Bazant (1994). Since the 
penetration load depends on the number of radial cracks, it is essential to understand 
how the number of cracks can be determined. In this dissertation, the crack initiation 
theory is stated in a form that facilitates further development. 

The structure of the dissertation is organized as follows. Chapter H related 
literature review. The basic concept and definition of crack initiation theory is 
introduced in chapter III, the problem of periodic parallel cracks initiated from the 
surface of a perfect brittle and elastic half plane under the action of initial strains will be 
discussed. The initial strains can be caused by a drop of temperature or by drying 
shrinkage of the material. The elastic stress distribution is solved by a integral 
equation. The solution of this problem can be related to the behavior of crack 
initiation that happens sufficiently far away from material boundaries. 

For the pavement problem, the cracking is normally confined in the top layer of 
the pavement, and the structure can be more appropriately modeled as an elastic beam. 
A new type of elastic analysis is introduced in chapter IV. Since the pavement layer is 
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usually about 10 inches thick, the nonlinear process zone of the crack is significantly 
large, and the consequence of the nonlinearity of the pavement material must be taken 
into consideration, which results in a further refinement of the crack initiation theory. 
However, the treatment should be considered as an empirical simplification. 

To fully consider the effect of nonlinear process zone, the cohesive crack model 
must be employed. In the cohesive crack model (CCM), the material does not fail 
immediately. Instead, there is gradual decreasing of bridging stress for increasing crack 
opening displacement, a phenomenon called strain softening, or simply softening. 
When the cohesive crack model is used, the concept of crack initiation needs to be 
generalized further. In addition, the theory of the cohesive crack model also needs to 
be expanded to consider the interaction between the cohesive cracks and the 
description of a unloading cohesive crack. Applying the cohesive crack model to study 
the problem of parallel cracks initiated from the surface of a half-plane is discussed in 
chapter V. 

Many problems of crack initiation theory still remain open. The main existing 
problems and possible solutions will be discussed. Future research direction will be 


described as a closure to this dissertation. 


CHAPTER n 


LITERATURE REVIEW 

In this chapter, some of the literature that are related to this study will be briefly 
reviewed. In addition, some nonlinear fracture models for quasibrittle material are also 
discussed because some of the analysis is based on these theories. To apply the crack 
initiation theory to concrete pavements, the literature of the concrete pavements design 
and research regarding temperature response are outlined. 

2.1 A System of Parallel Cracks from a Smooth Surface 

Consider the two-dimensional problem of a homogeneous isotropic elastic half- 
plane in which a system of parallel equidistant cracks normal to the surface is produced 
by cooling or drying shrinkage. This problem arises in many applications. In the later 
70's and early 80's an interest to extract heat from the hot-dry-rock by circulating the 
water through the cracks attracted a lot of researches. Also, there are applications of 
shrinkage cracking of concrete. Base and Murray (1982) called that "shrinkage 
cracking" is the most troublesome and frustrating quirk of structural concrete. 

Bazant and Ohtsubo (1977) first carried out stability analysis in which the 
conditions of stability of a system of Mode I cracks propagating along given paths 
were determined by analyzing the second variation of the work needed to create the 
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cracks as well as by formulating the conditions of adjacent equilibrium. The crack 
system that is in equilibrium state is stable if the second variation of the strain energy 
S 2 W is always positive, and unstable if it is negative, i.e. 

> 0 stable 

tfW = -jy£ t W (f Sa t Sa J < = 0 critical ’ (2-1) 

^ 1=1 ;=1 < 0 unstable 

Bazant, Ohtsubo, and Aoh (1979) continued their work to determine the critical 
states, bifurcation of equilibrium path and postcritical behavior. When a system of 
parallel equidistant cooling cracks propagates, it reaches a critical points, and the 
equilibrium path of the system bifurcates. The stable post-critical path consists of 
extension of every other cracks upon further cooling, initially with a crack jump at 
constant temperature, while the intermediate cracks stop growing and gradually 
diminish their stress intensity factor. Subsequently the leading cracks grow at equal 
length until they again reach a critical state, at which every other crack stops growing, 
and the process in which the crack spacing double is repeated. In both papers the finite 
element approach was adopted for completing numerical analysis. Bazant and Wahab 
(1979) also discussed the effect for different temperature profiles. 

Nemat-Nasser, Keer, and Parihar (1978) formulated the problem in terms of a 
singular integral equation of Cauchy type with continuous dislocation functions as basic 
unknowns for alternative cracks. They analyzed the crack growth regimes based on 
stability analysis of crack growth by checking the derivatives of the stress intensity 
factors with regard to crack lengths. It is equivalent to the Bazant's stability analysis. 
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They pointed out that the equal crack growth regime is stable as long as the following 
condition holds 


dK. dK 2 
— L = — - < 0 
da x da 2 


(2-2) 


Where K x and K 2 are stress intensity factors, a x and a 2 are crack lengths for the 
corresponding crack. An unstable critical state is reached when 


dK, SK 2 
= — 2-=0 

da x da 2 


(2-3) 


After this state, one crack, say crack 2, stops as the other crack, i.e. crack 1, grows 
spontaneously. 

Keer, Nemat-Nasser and Oranratnachai (1978) extended their work to the array 
whose unit cell contains three interacting cracks. The model can exam the crack 
growth regime at and after the second critical state, at which every other of those 
cracks which had continued growing stop growing, while the remaining one continue to 
grow at a faster rate. It is shown that when the temperature profile is in the form of an 
error function, the inclusion of the third interacting crack changes the previous obtained 
results qualitatively (i.e. no crack closure is attained in this case). The lack of 
symmetry in this model implies that crack extension involves in both Modes I and II, 
and the calculation requires the correspondingly great number of simultaneous integral 
equations, which causes the computational complexity. 


The above stability analysis requires accurate estimates of the values of the 
stress intensity factors at various crack tips, and their derivatives with respect to the 
crack lengths. Actually K\ is obtained by extrapolation of the solution of system 
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equations, and the calculation of 



is more difficult and less accurate. 


Sumi, 


Nemat-Nasser and Keer (1980) proposed a combined analytical and finite-element 
solution method which leads to a rather effective solution procedure. Nemat-Nasser, 
Sumi and Keer (1980) also used this method to study stable and unstable bifurcation 
points and critical points. In particular, they pointed out the stability of this kind is 
highly imperfection-sensitive. 

To verify the previous theoretical results Geyer and Nemat-Nasser (1982) 
conducted experimental investigation on thermal induced parallel edge cracks in half- 
plane of brittle material. Glass plates were heated to a uniform temperature and then 
brought in contact with dry ice. The thermal contraction of boundary layer produced 
interactive tension cracks. They observed that the cracking occurred in a dynamic 
fashion seconds after initial contact with the dry ice bath. However, dynamic crack 
growth is related more to the problem of ‘brack initiation” than ‘track growth”, as 
was recognized by the authors. 

Bazant and Wahab (1980) used previous formulation of stability condition to 
investigate the reinforced concrete (one layer of steel reinforcement) by finite element 
approach. The bond slip (between concrete and steel) length was assumed and taken 
into account. 

In Chapter 12 of his book (Bazant and Cedolin,1991), Bazant summarized the 
stability analysis of a system of parallel cracks. In all these works it was found that 
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instability of a system of parallel shrinkage or cooling cracks may cause some cracks 
to close and the remaining ones to extend and widen. 

However, most attention is focused on the problem of how cracks develop in 
structures that already have them. Since cracks propagation happens after crack 
initiation, thus the problem of crack initiation becomes an issue. 

Although there were some crude estimates. Keer, Nemat-Nasser (1978, 1979) 
estimated crack spacing with several approaches. One is by energy consideration, its 
idea is that the portion of the total strain energy before cracking is used to generate 
new surfaces as the thermal crack initially form ( GW = 2 ay). Then the minimum crack 
spacing b is obtained in this way as 

4 =*»i§ < 2 - 4 > 

where a is crack length, d is penetration depth of drying or cooling, and 0 < 6< 1, 
denotes the fraction of the total strain energy that will be released. The coefficient k 0 
depends on the temperature profile and material properties. However, neither ratio a/d 
nor fraction 6 can be determined in their analysis, as was admitted by Nemat-Nasser 
(1979). Also Bazant and Wahab (1979) and Bazant and Cedolin (1991) gave a similar 
estimate of the minimum cracking spacing for the parabolic temperature profile as: 

20(1 - v)y a 

0(1 + v)c?£s.T 2 E d K ) 

where y is the surface energy density, a the coefficient of thermal expansion, AT the 
surface temperature drop. 
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2.2 Classical Cohesive Crack Models 

Barenblatt (1959a, 1959b, 1962) was the first one to challenge the concept of 
stress singularity in fracture mechanics. His basic argument is that stress should always 
be finite. His mathematical model of crack differs from that of Griffith's model in that 
the cohesive interatomic and/or intermolecular forces are explicitly included in the 
equilibrium of the crack. It is noted that in Griffith's approach the cohesive forces are 
only accounted for in surface energy computation but ignored totally in the equilibrium 
consideration. Barenblatt pointed out an interesting fact that the smoothness of crack 
opening profile is equivalent to the finiteness of stress at crack tip, which can be seen 
from the solution of the infinite plane with a semi-infinite cut on the negative x-axis. 
Besides the external loads that make the crack open, there are distributed forces p(x) 
defined symmetrically on the both sides of the crack that tend to close the crack. The 
normal stress and displacement in y-direction near the crack tip can be expressed as 

K +p(0 ) + O(^); v _ T j ( il , ^ + 0( ^ ) (2 . 6) 

where s, is in the positive x-direction while Sj in the negative x-direction, the sign in 
front of displacement indicates that each crack face goes to opposite directions, v is the 
Poisson's ratio. When K is zero in these expression, then crack opening is smooth at 
the crack tip. Besides, it is noted that when stress singularity is nullified, the normal 
stress component in y-direction is continuous at crack tip. The connection between the 
smoothness of crack opening and the finiteness of stress was perhaps first demonstrated 
by Westergaard (1933) in a paper on the stress distribution in a RC beam with crack. 
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The cohesive force generally depends on the crack opening displacement, 
therefore the cohesive crack models are generally nonlinear. However, Barenblatt 
argued that if the process zone is small when compared with the crack length, then the 
cohesive force can be considered to be independent of crack opening displacement. 
Furthermore, the cohesive force can be used to determine the modulus of cohesion K,.: 


_1 fP(s 2 )ds 2 


K ‘ =1 J 

Po 


where lp is the length of process zone. It is to be remembered that now p(x) depends 


only on its position rather than crack separation. The distinction between the classical 


cohesive crack model and the cohesive crack model discussed herein lies in whether 


cohesive forces are considered dependent on the evolution of crack opening 
displacement. By adopting the above assumptions, the process of crack formation is no 
longer considered. 

It should be noted that in the original definition, the factor (l/7t) is not included. 
Furthermore Barenblatt argued that the modulus of cohesion has the same relation with 
the strain energy rate G as the stress intensity factor K when crack is about to 
propagate. It is therefore fair to say that Barenblatt's theory did not yield any 
operational difference from linear elastic fracture except for a more realistic picture of 
what is happened near crack tip. 

The postulation of equivalence between Barenblatt's criterion and Griffith's 
criterion was examined mathematically by Willis (1967). Willis found that in order for 
the two theories to be compatible, it is necessary to require that the cohesion zone be 
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small compared to the crack length, and the load level must be very low with respect 
to the Young's modulus of the material. The problem of the profile of crack separation 
within cohesion zone was addressed by Smith (1974). Instead of solving for crack 
separation for given cohesion law, which is a nonlinear problem, an inverse approach is 
employed, that is, stress-separation law is deduced from given separation which is 
usually of very simple form. At this point, one should keep in mind that Smith's results 
rely completely on the assumption of smallness of cohesion zone, any extension to 
nonlinear fracture mechanics (e.g. Karihaloo et. al. 1989) is, therefore, unlikely to be 
justifiable. 

Dugdale (1960) was credited for demonstrating the usefulness of cohesive 
crack model in nonlinear problem. Actually his objective was to find the extent of 
yielding in front of crack tip as a function of loading over yielding stress. Assuming that 
the plastic region is a thin strip extending from the crack tip, the cohesive force is con- 
sidered to be constant at the magnitude of yielding stress. Since no stress should exceed 
yielding stress, the stress intensity factor should be zero which in turn becomes a 
condition for determination of cohesive zone length. According to Dugdale's 
experiment, the theory provided excellent prediction. Eight years later Cotterell (1968) 
affirmed that Dugdale's model is equally applicable to organic glass. 

The drawback of Dugdale's work is that no crack propagation criterion was 
provided, which was not the intention of his paper, the cohesion zone can be stretched 
to arbitrary length for sufficient load level. To remedy the situation, the crack opening 
displacement (COD) was proposed to serve as the cracking criterion for Dugdale 
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model. When COD exceeds some critical value, which is assumed as a material con- 
stant, the crack will propagate. Later it will be seen that COD is not a constant in the 
cohesive crack model. 

2.3 Fracture Mechanics of Concrete 

Quasibrittle materials such as concrete, rock, ceramic, paper, lumber etc., are 
very different from metals, because they exhibit significant softening behavior while 
plastic hardening deformation is negligible. On the other hand, these materials are not 
as totally brittle as materials such as glass and mica. Quasibrittle materials also can 
include advanced cementitious composite, fiber reinforced composite, etc. (Shah 1991). 
In quasibrittle materials fracture is preceded by a fracture process zone, which may be 
caused by microcracking or other material defeats that the major source of nonlinearity 
come from. Some experiment observations of concrete indicate that within this zone, 
fracture is manifested in the form of aggregate debonding and overlapping 
microcracking, mainly parallel to each other, with some intact grain bridges between 
them (Mier, Rots and Bakker, 1991). Concrete has an added complication in that the 
scale of its mesostructure is relative large compared to other materials. Kaplan (1961) 
applied linear fracture mechanics to concrete for the first time. Since then most efforts 
have been directed towards better understanding of the pronounced nonlinearity in 
concrete and how to deal with such nonlinearity as has been summarized in the 
excellent review articles by Mindess (1983, 1984). 
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Kaplan (1961) was the first one to acknowledge the slow crack growth in 
concrete, which is a stable crack extension under increasing load before final rupture 
occurs. Kaplan was able to relate the length of slow crack growth with the size of the 
specimen. As a matter of fact, the phenomenon of slow crack extension is not unique to 
concrete. Irwin (1958) was probably the first one to undertake this problem by 
identifying it with plastic deformation. Slow crack led its way to the so-called R-curve 
approach by Broek (1968) in his study of an aluminum-copper alloy. Since then the R- 
curve approach became a very popular practice also in concrete fracture mechanics 
research (e.g. Bazant, 1984b; Wecharatana and Shah; 1983a, 1983b; Mai, 1984; Foote, 
1986, etc.). Although Broek (1986) realized that R-curve is not likely to be a material 
property, as has been indicated by numerous researches that R-curve is dependent on 
the loading configurations (Wecharatana, 1983b, e.g.) as well as specimen size ( Foote, 
et.al, 1986). However, due to the fact that there has been no better approach available, 
R-curve method still enjoys wide applications among some of more pragmatic 
researchers (e.g. Mobasher,1989; Shah, 1990). 

Another difficulty with linear fracture mechanics of concrete is the so called 
notch sensitivity problem. Griffith (1924) is the first one to consider the smallest size of 
crack which can cause the observable decrease in the strength of materials. For glass, 
Griffith found that this size is about 1.5*1 0* 3 (mm). When this size is relatively large, 
the material may be called notch insensitive. While many researchers agree that 
hardened cement paste is a notch sensitive material, the opinions about concrete and 
fiber reinforced concrete are divided, both sides claimed having experimental evidence 
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backups. When it comes to the experimental measurement of K,. and G c , the results 
are again highly contradictory, especially in the study of the dependency of and G c 
on the geometry of specimens. Although obviously the highly heterogeneous nature of 
concrete is responsible at least in part for this apparent inconsistency, nevertheless by 
the end of seventies the linear fracture mechanics was not considered directly applicable 
to concrete material by most researchers. 

Hillerborg's fictitious crack model (FCM) was advanced during the late 
seventies and early eighties (1976,1985) as the first nonlinear fracture mechanics model 
in concrete. Apparently encouraged by the success of Dugdale's model, Hillerborg 
proposed to use a declining cohesion law (i.e. stress-separation relation) to describe the 
strain softening behavior of cementitious materials. The total area below the stress- 
separation curve defines a quantity called fracture energy G fi which is experimentally 
found to be less dependent on specimen size than any other fracture parameters 
(Hilsdorf and Brameshuber, 1984). An important difference from that of Dugdale's 
model is that the cohesive stress now is displacement (crack separation) dependent, 
therefore the model became truly nonlinear because the cohesive stress is neither 
constant nor proportional to the crack separation. It is so difficult that no analytical 
solution even for the simplest geometry has been found. Convenience in finite element 
application is one of the objective in the original paper (Hillerborg, 1976), since the 
removal of stress singularity made it possible to use regular coarse mesh. 

When Bazant proposed his crack band theory, his major objective seemed to 
aim at the formulation of a nonlinear fracture model for concrete (and its similar 
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material, of course) such that the distributed microcracking as well as coalesced (that 
is, localized) cracking can be treated in a unified way (e.g. 1984a, 1984c, 1986a, 
1990a, 1990b). The starting point is a numerical technique called smeared crack 
simulation proposed by Rashid (1968) where crack is represented by elements with 
tensile Young's modulus equals to zero. Bazant employed the method to include strain 
softening in the formulation. No sooner than the scheme was implemented, the result 
was found mesh dependent, i.e., the energy dissipated by strain softening diminishes as 
the mesh shrinks. The strategy used in his crack band model is to take softening 
modulus E t as a function of Young's modulus E, tensile strength /„ the fracture energy 
Gy- and the cracking front width w e such that the fracture energy Gy is constant: 

G f = (^ +-^~) w c = constant (2-8) 

Z fc, -Jb T 

where the factor that multiplied by w c can be recognized as the specific strain energy, 
that is, the area under stress-strain diagram. But if E r is determined independently, then 
w c can be determined from the above equation. Bazant considered w c as a material 
property and refused to take elements smaller than w c . 

A comparison between FCM and crack band model was made by Bazant 
(1986a) in which Bazant recognized certain similarity after microcracks are localized 
into one macroscopic crack. It is fair to say that Bazant's goal is more ambitious, it is a 
serious strive towards a unified theory of distributed cracking as well as localized 
cracking. On the other hand, Hillerborg's objective is modest, he chose to deal only 
with those cases where localization is known a priori. Consequently FCM is simple 
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both mathematically and physically. The introduction of stress-separation relations in 
place of stress-strain relation circumvented the predicament of unobjective softening 
modulus E t , leading to a much simpler numerical scheme. The attractiveness of this 
approach was enhanced by Ottosen (1986). Using stability and uniqueness arguments, 
he demonstrated the ominous incongruity accompanied with strain softening concept. 

One of the interesting feature of FCM is the unified treatment of notched 
fracture and unnotched fracture. As it was mentioned in the beginning of this section, 
notch is not necessary for concrete to rupture by cracking as most metals do, since it is 
its only way of rupture under tension. Extensive calculations of FCM on 3-point beams 
was made by Carpinteri (e.g. 1989) for different sizes and different notch ratios 
including zero notch ratio. A boundary element simulation was made by Liang and Li 
(1991a). The size effect of FCM is examined (Liang and Li, 1991b). It is fascinating to 
see how size effect curves gradually shift towards a curve of bending strengthening. Li, 
et.al. (1986) used the Green's function method to solve FCM for a centrally cracked 
infinite plane using different stress-separation curves with the same G f and the same 
tensile strength. Among many conclusions, they found that the shape of the stress- 
separation curve has a significant effect on the load versus crack opening diagram, and 
that the length of process zone is not a material property, but rather a variable 
depending on the loading configuration and structural geometry. Liaw et.al. (1990) 
recently demonstrated that better results can be obtained by adopting a more refined 
stress-separation curve in FCM computation. Li and Liang (1993) developed a theory 
of CCM in which the peak load of the Griffith problem (a Central tensile crack in an 
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infinite plane) was solved through the condition of stability limit. The critical 
condition can be transformed into a linear eigenvalue problem under the assumption of 
linear softening law. Using the same technique, Li and Hong (1992) solved cohesive 
crack problems of double-notched or center notched infinite strip under remote tensile 
loading. 

There is another nonlinear fracture mechanics model of concrete material called 
two-parameter model by Jenq and Shah (1985a,1985b). As is well known that the 
critical stress intensity factor K c and critical energy release rate G c calculated from the 
peak load depend on the dimensions of specimens, Jenq and Shah proposed to use 
crack mouth opening displacement (CMOD) as a second crack parameter to 
complement linear fracture mechanics parameter such as K c . In computation, CMOD is 
used to determine the length of process zone in which a stress-separation relation is 
assumed. After the adjustment caused by the bridging effect in the process zone, the 
calculated net stress intensity factor at peak load is found to be constant by their 
experiments. Actually the idea of using cohesive force as a modifying factor in linear 
fracture mechanics exits in Shah's research before 1985, as can be seen in their earlier 
papers by Wecharatana and Shah (1983 a, 1983b), and a discussion by Hillerborg 
(1984) in which some ambiguity in definition as well as in concept was commented. 
The two-parameter model published later is just a manifestation that Shah hardened his 
position in defending his notion that the stress singularity still exists in spite of the 
presence of pronounced cohesive force in process zone. Recently Shah in a keynote 
lecture (1989) vindicated his model by categorizing the dissipation mechanism in 
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concrete into three sources, one that corresponding to surface energy, one due to 
microcracking and the one that dissipated in the wake of crack path. According to this 
classification Shah claimed that the two-parameter model is the best in terms of energy 
balance consideration. 

It is a generally accepted idea that when concrete specimen is small, the 
strength criterion will dictate the rupture load, and when the dimension becomes 
sufficiently large, the well-known inverse-square-root scaling law of linear fracture me- 
chanics will govern. Now the question is what happens in between these two extremes. 
Walsh (1972) tried for the first time to establish the size effect. It is very interesting to 
note that he used (f/aj 2 as ordinate and d/d a as abscissa to plot his result which is a 
straight line, and this is basically the same form as the one proposed by Bazant (1984a) 
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2.4 Crack Spacing in Concrete Pavement Due to Temperature Effects 

The first Portland cement concrete pavement in the United States consisted of a 
10-ft wide by 220-ft long slab constructed in Bellefontaine, Ohio in 1891. During the 
following 100 years tremendous increased pavements were build at various locations in 
the Unite States. In 1991, the Federal-aid highway system comprised more than 
850,000 miles of pavements of various design. (FHWA, NCHRP Synthesis 189, 
1993). 

It is necessary to briefly review the basic consideration regarding temperature 
effect in concrete pavement design from “AASHTO Guide for Design of Pavement 
Structures 1993”. The rigid pavement (Portland cement concrete pavements) includes 
plain jointed (JCP), jointed reinforced (JRCP), and continuously reinforced (CRCP). 
Joints are placed in concrete pavements to permit expansion and contraction of the 
pavement, thereby relieving stresses due to environmental changes (i.e., temperature 
and moisture) and friction. There are three types of joints for JCP: contraction, 
expansion and construction and their functions are as follows: 

(1) Contraction or weakened-plane joints are provided to relieve the tensile stresses 
due to temperature, moisture, and friction, thereby controlling cracking. If contraction 
joints were not installed, random cracking would occur on the surface of pavement. 

(2) The primary function of an expansion joint is to provide space for the expansion 
of the pavement, thereby preventing the development of compressive stresses, which 
can cause the pavement to buckle. 

(3) . Construction joints are required to facilitate construction. 
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For contraction joints, according to the AASHTO (1993), the spacing 
decreases as the thermal coefficient, temperature change, or subbase frictional 
resistance increases; and the spacing increases as the concrete tensile strength increases. 
The spacing also is related to the slab thickness and the joint sealant capabilities. It is 
suggested, as a rule of thumb, the joint spacing (in feet) for plain concrete pavements 
should not greatly exceed twice the slab thickness (in inches). For example, the 
maximum joint spacing for an 8-inch slab is 16 feet. In ‘Rigid Pavement analysis and 
Design’l(FHWA-RD-88-068, 1989) joint spacing for plain concrete pavement ranges 
between 12 and 20 feet, with maximum 20 feet. California instituted the practice of 
specifying joints at 12, 15, 13, and 14 feet. The use of a random spacing pattern was 
instituted in the early 1960’s to minimize excessive vibration problems (FHWA-RD-86- 
040, 1986). 

The width of the joint is controlled by the joint sealant extension. The depth of 
contraction joints should be adequate enough to ensure that cracking occurs at the 
desired location rather than in a random pattern. Normally, the depth of transverse 
contraction joints should be 1/4 of the slab thickness. These joints may be developed 
by sawing, inserts, or forming. Time of sawing is critical to prevent uncontrolled 
cracking. 

The use of expansion joints is generally minimized on a project due to cost, 
complexity, and performance problems. They are used at structures where pavement 
types change (e.g., CRCP to jointed), with prestressed pavements and at intersections. 
In ‘Pavement Structural Design Practices” (1993) it is stated that the use of expansion 
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joints was virtually standard practice by the mid 1930s, but has since been found to 
be unnecessary except where the pavement abuts a structure. 

Major distress problems of concrete pavements generally start with crack 
formation caused by the combined effects of traffic load and service temperature. 
Water and salt can easily infiltrate into the pavement at the location of cracks and 
create durability and structural problems. 

The structural response of concrete pavements under traffic loads is highly 
dependent on temperature and its variation. Temperature effects on concrete pavement 
behavior have been recognized since the mid-1920s. Westergaard (1926) identified 
temperature curling as an important parameter affecting the structural behavior of 
concrete pavements. Westergaard’s method (1925) for computing stresses in concrete 
pavement was based on an assumption that the subgrade acts as a Winkler foundation: 
the pressure between slab and subgrade is proportional to deflection. Teller and 
Sutherland (1935) reported the results of tests conducted on concrete pavements to 
study the effects of variations in temperature and moisture. Lang (1940) studied the 
movement of concrete pavement slabs resulting from changes in temperature and 
moisture. Friberg (1954) presented a mathematical evaluation of horizontal slab 
movements, and effect of the subgrade frictional resistance, on stress development in 
long pavement slab. Harr and Leonards (1959) conducted laboratory tests to measure 
temperature curling and compute subsequent stresses. They correlated the results from 
the laboratory tests with predicted response with an analytical model that they 
developed. 
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Armaghani, Larsen and Smith (1987) try to more precisely describe the 
displacements (vertical and horizontal) of a concrete pavement slab associated with 
temperature variation and weather. Temperature data, accumulated by Bergen from 
1983 to June 1986 from a test road (in Florida), are analyzed. Slab displacements were 
monitored and evaluated. It is found that maximum daily displacements were 
concurrent with maximum temperature differentials in the slab. They pointed out ‘in 
many analysis of thermal stresses, temperature gradients have been assumed to be 
linear. This assumption has simplified the modeling of pavements without significantly 
affecting the accuracy of the computations. Therefore for all practical purposes the 
temperature gradient can be approximated by a linear curve.” Richardson and 
Armaghani employed a parabolic function to model the nonlinear temperature parallel 
to the pavement thickness. In their particular case the nonlinear temperature stress is 
only 17 percent of the flexural strength. Although the importance of nonlinear 
temperature may not be important for deformation due to normal daily temperature 
fluctuation, they did not rule out the importance of the nonlinear temperature 
distribution during dramatic temperature changes. Because cracks are most likely to 
initiate during dramatic temperature changes, the nonlinear distribution can be 
important for our purposes. 

In 1970’s, continuously reinforced concrete pavement (CRCP) increase 
dramatically in use, for instance, from 1961 to 1971, increase is 20 time. In design of 
CRCP, allows the effects of shrinkage and temperature change to produce random 
cracks but keep the cracks tightly closed together. Comparing with plain concrete 
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pavements, it has much shorter crack spacing and width. The average mature crack 
spacing is around 5 feet. The average percentage of steel is 0.6 (the more steel, the 
shorter the crack spacing.) 

Most studies on the crack formation mechanisms in concrete pavements are 
still limited to the conventional models developed on the basis of stress-based or strain- 
based elasticity analysis. Jenq, Liaw and Kim (1993) applied the fracture mechanics to 
study the effects of temperature on early crack formation. To properly control the 
occurrence of random cracking, saw-cut grooves are generally introduced at the 
earliest possible age of the concrete pavement. Adequate groove depth must be 
provided to ensure that the transverse cracks will be confined at the location of the 
groove. The cohesive crack model and finite element method are used, once knowing 
the expected temperature deferential and thickness of the pavement, they can determine 
the timing and groove depth. Based on their theory the spacing of saw-cut groove can 
not be determined. 

It is reported by Federal Highway (FHWA, 1990) that thermal cracking of 
asphalt pavement continues to be a problem in many parts of the Unite States. The low- 
temperature shrinkage cracking is concerned in Canada and northern United States, 

and thermal fatigue cracking is now recognized as a problem in more temperate 
climates. The early models to predict thermal cracking are based upon empirical or 
statistical relationships that relate cracking to various asphalt specification data and 
environmental parameters. 
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More recently, fracture mechanics theory was used to develop a computer- 
based model that can be applied to the thermal cracking problem. During the 
development of this model, it was necessary to resort to statistical regression equations 
to predict fracture properties, rather than to incorporate fundamental asphalt or mixture 
fracture properties into the model. Thus, neither the statistical based models nor the 
more recently developed fracture asphalt model is a completely mechanistic approach 
for relating fundamental asphalt or mixture properties to the incidence of thermal 
cracking. 

The study of this dissertation treat asphalt concrete as a quasibrittle material 
with significant post-peak softening behavior. Actually asphalt concrete is more like a 
viscoelastic material in the warm temperature. However, during winter weather, 
viscosity of the asphalt concrete become less important, and as an approximation it can 
be reasonably excluded in our analysis. Although the materiel properties are 
temperature dependent, the dependence is not as pronounced if temperature is low. 

The cracking of the asphalt pavement is sometimes measured by the number of 
cracks per mile. Some data are reported (Ruth, 1982) as to how this cracking index 
increases year by year for different type of asphalt concrete mix. However, in this 
report, only the cracks that cross the whole traffic lane are counted. Conceivably, if 
shorter cracks are also included, the value of the crack index should be multiplied by a 
factor of 2 - 4. Accordingly, the first year average crack spacing is about 30 meters if 
only full grown cracks are counted. If shorter cracks are also included, then the 
spacing is about 7 to 1 5 meters. 


CHAPTER m 


INITIATION THEORY — FOR BRITTLE MATERIALS 

Because of cooling, externally applied loads, residual-stress build-up due to 
creep, loss of moisture and consequent shrinkage, or other natural of imposed 
processes, cracks often form from smooth surfaces. Problems of this kind include 
shrinkage cracks in drying concrete, Shrinkage cracks in polymers due to aging and 
loss of moisture, desiccation cracks in deserts and at the bottom of dried up lakes, to 
name just a few. The situation may be idealized as a system of parallel Mode I 
equaldistance cracks normal to the smooth surface of half-plane. 

The crack propagation of a system of parallel cracks was studied in detail with 
respect to a proposed hot-dry-rock geothermal energy scheme in late ‘70s. The stability 
analysis of crack propagation was done by Bazant and Ohtsubo (1977), Bazant, 
Ohtsubo, and Aoh (1979) and Bazant and Wahab (1979) that they adopted the finite 
element approach; also by Nemat-Nasser, Keer, and Parihar (1978), Keer, Nemat- 
Nasser, and Oranratnachai (1979), Sumi, Nemat-Nasser, and Keer (1980) that they 
formulated the problem in terms of a singular integral equation. But most attention is 
focused on the problem of how cracks develop in structures that already have them 
(crack propagation happen after crack initiation). Although Keer, Nemat-Nasser 
(1978, 1979) gave some crack spacing estimate based on arbitrary assumption. 
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The purpose of this research is to provide a theoretical analysis and numerical 
approach for the crack initiation problem. The initiation of thermally induced equally 
spaced and equal lengths parallel cracks in a half-plane consisting of a homogeneous 
isotropic linearly elastic brittle material are studied. The problem is stated and basic 
model is outlined in section 1. In section 2 the basic three conditions for governing 
crack initiation of brittle material are cast. The corresponding thermoelasticity problem 
is formulated in the form of integral equations in section 3, and the numerical method 
of solution is outlined in section 4. Numerical results, relevant discussion and 
comparison with other theoretical results at limit cases are then presented in section 5. 
The experiment evident of glass comparing with the solutions are reported in section 6. 
Some comments and a few conclusions are presented in section 7. The possible 
experiment is suggested in the final section. 

3.1 Statement of Problem 

Consider an elastic half-plane of unit thickness that occupies the region y > 0 
(Fig.3-la) and is under plane stress. The material is homogeneous and isotropic. The 
half-plane is subjected to initial strain e x = e 0 f (y/d) which may be caused by cooling 
or drying of the surface; / is a given strain function and d is the penetration depth of 
cooling or drying. The initial stress is o x = Ee x = SgEfly/d). The general shape of 
function fly Id) is shown in Fig.3-lb. The initial strain profile may be considered either 
an error function (which is the exact solution of the linear diffusion problem) or a 
parabolic function (which is an often used approximation). These functions are 
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( 3 - 1 ) 


/O0 = i — j= \f ye ~“ 1(/u or / 0 ') = 0-> / ) 2 o <3/ < i 

with maximum value /(O) = 1. We shall use both functions. Basically they give similar 
results. 

The phenomena of crack formation are often complicated by fhe randomness of 
material inhomogeneity. The initial cracks rarely appear in a regular and systematic 
way. In the case of large concrete blocks, the cracks are never straight and the spacing 
between these cracks is hardly uniform. In the case of a drying lake bed, the crack 
pattern may not be exactly hexagon; it may be in the shape of a pentagon or heptagon. 
However, it is also true that hexagon is the most common shape, and the sizes of these 
hexagons are almost uniform. In the case of penetration of a sea ice plate, it is found 
that the total number of radial cracks varies from experiment to experiment, but 
generally larger punch sizes produce more radial cracks (Frankenstein, 1963). Despite 
the randomness, there must nevertheless be deterministic laws underlying these 
phenomena. 

However, due to lack of information, and because of the accompanying 
analytical difficulty, we dramatically simplify the problem into the idealized case of 
cooling of a perfectly homogeneous and isotropic elastic half-plane, in which we 
expect the initial cracks to be straight, perpendicular to the free surface (system of 
parallel cracks), their spacing uniform and their pattern regular. The temperature profile 
in the solid is not altered by the formation and extension of these cracks. Moreover, to 
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simplify the analysis further, we assume that the initial cracks are of equal length a in 
this chapter. Thus we may concentrate on one cell of width b as shown in Fig.3-lc. 

The problem is first decomposed into two separate problems. The first is the 
structure without cracks with the'initial (or residual) stresses applied. The second is 
the structure without initial stress, but with cracks and the crack surface traction that 
are equal in magnitude to the initial stresses but opposite in direction. The solution is 
obtained by the superposition of these two problems. Since there is no stress intensity 
factor in the first problem, we can, therefore, study only the second problem. This is the 
standard approach studying the fracture mechanics problem in elastic materials with 
residual stresses. 

3.2. Conditions Governing Crack Initiation from Smooth Surfaces 

It is well known from experience that the formation of initial macroscopic 
cracks from a smooth surface of brittle material is a sudden event. The phenomenon of 
crack initiation involves transition between two states: One is the initial equilibrium 
state in which there are no cracks. The other is the first state of stable equilibrium at 
which the cracks have formed and have come to a stop but are simultaneously in a 
critical state from which they can propagate further in a stable equilibrium manner. In 
between these two states, the cracks are unstable and grow dynamically. The reason 
that the transition is dynamic is that the energy release rate as a function of crack length 
begins increasingly from zero. At the same time, since we treat the problem according 
to linear elastic fracture mechanics (supplemented by the strength criterion), the 
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fracture energy is constant and nonzero. This is what causes the initial instability of 
crack growth. In this study we are interested only in the two equilibrium states and the 
basic laws governing them. Because the material is considered to be elastic the 
response is path independent, and so the dynamic transition between these two states 
need not be analyzed in order to solve the crack spacing problem. 

The first condition governing the crack initiation is a stress condition: The 
tensile stress caused by load before cracking must reach the tensile strength of the 
material at least at one point. For any load below this level the material simply can not 
start to break. 

The second condition is that the final state, as we have defined it, satisfy the 
Griffith law that the energy release rate be equal to its critical value. If the energy 
release rate were above the critical value, then the cracks already formed would be 
unstable and would not stop at the final state. On the other hand, if the energy release 
rate were below the critical value, then the cracks would have had to stop earlier and 
that state would not represent the final state. 

The third condition is provided by the law of energy conservation, which 
requires the potential energy of the structure released due to crack jump to be fully 
converted into surface energy of the newly formed cracks. Note that the Griffith law 
(the second condition) is also a statement of energy conservation which is however 
applied to an infinitesimal crack extension and is represented by a differentiation of the 
potential energy, whereas the third condition is energy conservation for a finite crack 
jump and is represented by a finite difference in potential energy. 
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The three conditions of crack initiation can simply be written as 

C7 x (x,0)>/, ; G(a,d,b) = G f ; U(0,d-,b)-U(a,d-,b) = aG f (3-2) 

where f( = tensile strength of the material; G/= fracture energy of the material defined 
as the energy required to open a unit crack fully, U = strain energy within the cell of 
width b. For simplicity, we use here for the stress condition the equality a x = ft, and 
postpone the consideration of the case of inequality. In the following derivation, we 
will make use of Irwin's formula for energy release rate in terms of the stress intensity 
factor, 


^ dU K 2 


(3-3) 


where K = mode I stress intensity factor at the crack tip , defined as 
K = lim a x [2 7t(y - a) ] 1/2 ,and E = Young's modulus of the elastic material. With this 
relation and based on the third condition, we can express the strain energy change as 

0 0,d;b)-U(a,d;b) = ±\‘K 2 da (3-4) 

E J0 


Furthermore, by combining the second and third conditions and rearranging, we can 
cast the three basic conditions as: 


cr x =/ f ; K 2 = EG f ; aK 2 =\°K 2 da (3-5) 

Since all the equations refer only to a unit cell of width b (crack spacing) the second 
condition of crack initiation furnishes a relation between the ratios of a/b and d/b. In 
this way, solution of the crack initiation problem yields K as a function of a/b and d/b 
for the given initial strain profile. In other words, there is a one-to-one correspondence 
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between crack spacing b and the loading depth d. Once the penetration depth d is 
given, the crack spacing b and the initial crack length a can also be determined. We 
will come back to this point later. 

3.3 Mathematical Formulation 

Since all these quantities are complicated functions of the geometrical 
configuration and the loading profile, the solution has to be numerical. 

The problem of a half space weakened by equidistant parallel surface cracks has 
been studied by Bazant et al. (1977, 1979) using the finite element method, by Nemat- 
Nasser, Keer and Parihar (1978) using the dislocation representation, and by Nied 
(1987) using the displacement jump as the basic unknown. The problem was also 
reviewed by Bazant and Cedolin (1991). The stresses for one normal surface crack in a 
half space were given by Keer and Chantaramungkom (1975), and Nemat-Nasser, Keer 
and Parihar (1978) modified the expression to obtain the stresses for an array of 
equidistant cracks on the surface. When all cracks are of equal length, as assumed 
here, the equation can be simplified. The condition that stress at the surface reaches 
the tensile strength may be written as 

<r,(0,y)=j-£m i,g«,»b,y)d< = -/, f(y/d) ( 3 -«) 

n =- oo 

where E is the Young's modulus for plane stress condition and changed to E / (1 - v 2 ) 
for plane strain condition, where v is Poisson's ratio. Function D(t) represent the 
dislocation density . The kernel function g can be expressed as 
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2(y + 0 (y + 3Q[(y + 0 2 -x 2 ] t 4 ty(y + 1 )[(y + 1 y - 3x 2 ] 

(y + 1 ) 2 + x 2 [(>> + 1) 2 + x 2 f [(y + 1) 2 + x 2 ] 3 

(y-t)[(y-t) 2 + 3x 2 ) 

[(y-t) 2 +x 2 ] 2 

(3-7) 


Since g is an even function of x, the series in (3-6) may be written as • 

f^gO,nb,y) = g(t,0,y) + 2^g(t,nb,y) (3-8) 

rr=-ao n=l 

Using the following two formulae (Gradshteyn and Ryzhik, 1965, pp. 23, pp. 36), 

(3-9) 
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the infinite series can be summed into a finite expression: 
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Also from (3-7) we have 
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It is further convenient to represent the kernel function in the following form 


g(t,0,y) + 2£g(f,wZ>,y) = f g, 
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where the function g { can be written as 
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Because function D(i) is singular at the crack tip, it is now convenient to introduce a 
new unknown function C(t) : 


D(t) = 
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(3-15) 


which is a smooth function. Equation (3-6) is thus transformed into the following 
equivalent integral equation for the unknown function C(t)\ 
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Also the unknown function C is a bounded function. Once the unknown 
function is solved, the stress intensity factor can be calculated as 
K= \m^27r{y-a) cr r (0 ,y) 
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(3-17) 


The detail proof can be found in appendix 1 . 


3.4 Numerical Method 

Since the first term of the expression of function g x is l/(t-y), (3-16) is actually 
a Cauchy's singular integral equation. To solve it numerically, we first normalize the 
interval (0,a) by defining 


s = t/a , x = yja, C(as) = f?(s) 

Now the singular integral equation (3-16) can be written in the form 

E f . = 
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(3-18) 


4 ^^°(l - 5 2 )' /2 * 


0 < x < 1 (3-19) 


Next, (3-16) is extended into the interval (-1, 1) by an even continuation: 

B(-s) = B(s) 0<s< 

In this way, (3-16) can be equivalently expressed as 
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With the Gauss-Chebyshev quadrature and the collocation technique as described by 
Erdogan, Gupta and Cook (1972, pp. 380-381, where n is replaced here by 2n+l), we 
can convert the singular integral equation into a discrete form as 
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7 = 1,2 ,n 


where n is the number of integration points and 
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v> yj = 
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then equation (3-22) may be written as 
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7=u " (3 - 25) 


where 


H(t,y, e) = ne{4t + 2 coth(y + t)-(y + 3t) coth 2 (7 + 1) 
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If we introduce 
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Then (3-25) become 
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Such a system of linear equations can be easily solved, for example, by the method of 
triangular factorization. Once we know A, we can calculate the stress intensity factor 
K in the form of dimensionless stress intensity factor N as 
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( 3 - 29 ) 


Note that the numerical solution does not yield directly A( 1). The closest data point on 
which the unknown function is defined is = cos(7r/(4/tf-2)). In theory, an 
extrapolation, for instance a quadratic extrapolation, is needed to find A( 1). However, 
when n is large enough (for instance, n — 50), the difference between ,4(1) and A{s , ) is 
in the third of fourth digit, and thus is negligible. 


The stress intensity factor, either in the form of K or the dimensionless form N, 
is obviously a function of a/b and d/b as well as the initial strain profile f(y/d). Our 
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purpose in the next section will be to find the relation between these geometric 
characteristics. 

In terms of dimensionless stress intensity factor N f we can rewrite the last 
equation in (3-5) as 


£ N\x, A )dx = eN 2 ( e , A) (3-30) 

Note that this equation is purely geometric, that is independent of the material 
properties. From this equation, we can establish a relationship, which turns out to be 
one-to-one, between e and A . Furthermore, let 


/„ = 


G f E 

IF 


or- 


f? 


(3-31) 


This frequently used material characteristic, which will be called the effective length (or 
characteristic length), is an important quantity in this analysis. (For the concrete with 
Young’s modulus E = 4.2 x 10 6 psi, fracture Energy Gf = 85 N/m = 0.485 lb/in, and the 
tensile strength f t = 3.45 Mpa = 500 psi, the effective length l 0 would be 8.15 inches.) 
All the lengths characterizing the geometry will be normalized with respect to l 0 . The 
second equation of (3-5) can now be converted to the form: 


N 2 {e,A) = Fo_ = —L^ 
v 2 nb 2 nb 


(3-32) 


where b * = b/l 0 is the dimensionless crack spacing. Also, the dimensionless crack 
length and loading depth can be defined as 


* ^ ^ L * 1 * 

=n b = eb 


d' = - = -b" = Ab* 
lo b 


a 


(3-33) 


39 

In the calculations, the ratio A is given as an input. Then, using (3-30) we can solve 
for e. After e and A are known, the value of b* is determined from (3-32). The 
dimensionless quantities a* and d* can be determined using (3-33). In this way, the 
problem is solved accurately and efficiently. 

The detailed procedure are given step by step as : 

(1) . give ratio A (penetration depth d over crack spacing b ); 

(2) . select ratio e (crack length a over crack spacing b)\ 

(3) . solve for ,4(1) through eq.( 3- 28) and calculated by eq. (3-29); 

(4) . check equation (3-30). If N 1 satisfy it go to next step otherwise go back 

to step 2; 

(5) . determine b*, a* and d* by using eqs. (3-32), (3-33). 

3.5 Analysis of Numerical Results 

Fig. 3 -2 shows the relation between a* and d*, as well as between b* and d*. It 
should be emphasized that the relation plotted in Fig. 3-2 between the penetration depth 
d* and the crack spacing b* (as well as the initial crack length a *) does not apply for 
the subsequent crack evolution. Rather, each point in Fig.3-2 represents an event of 
crack formation. After the crack is formed, there are other laws that govern the further 
growth of the crack system, which have been discussed in detail by Bazant et al. (1977, 
1979, 1991) and Nemat-Nasser et al. (1978, 1979). 

The solid curve shows the error function and the dashed curve the parabolic 
function. As can be seen from Fig.3-2, the difference in the final results between these 
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two profiles is not significant. Thus, all the following analysis refers to the parabolic 
profile only. As the dimensionless loading depth d* increases, the spacing b*, starting 
from a very large value (which is actually infinite), decreases monotonically. However, 
the initial crack length first decreases. After it reaches its minimum value, it increases 
with d* towards infinity. Such a behavior must reflect the requirement of energy 
balance. The crack driving force is controlled by the loading depth d*. For small d* 
there is only a small amount of energy available, and so the cracks must be very sparse, 
crack spacing b * must be very large and the crack length a* must be very small. On 
the other hand, a larger d* provides a larger amount of energy, and therefore a smaller 
crack spacing. Because the rate of decrease of b* is initially dramatic, the energy 
available for each crack must be reduced. That is why initially the crack length must 
decrease with the loading depth. After the rate of decrease of b* becomes less 
dramatic, the energy availability for each crack somehow catches up, and then the crack 
length a* begins to increase monotonically with d*. 

For very large d*, the crack length a* is also very large. Therefore, for the 
crack tip, the free surface at y = 0 is no longer important. In addition, when the loading 
depth increases unboundedly, the initial strain distribution becomes uniform. The 
problem is thus transformed into an array of semi-infinite cracks under the uniform 
surface pressure //. The stress intensity factor for this problem can be solved 
analytically using Fourier transformation method (Tada, 1985), and the result is: 




(3-34) 
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Substituting this solution into the second equation of (3-5), we obtain the theoretical 
result for the lower limit of crack spacing b* = 2. Such a limit serves as a check for our 
numerical calculation. Specifying a large value of d/b, we can solve for b* y which is 
found to be always larger, but very close to 2. For example, when d* — 10 4 , b* 
=2.026. If d* = 10 5 , then b* = 2.013. For the concrete with the effective length l 0 = 
20cm = 7.8 inches, the lower limit of crack spacing would be 15.6 inches. 

The other limit corresponds to a small d* value when the spacing is infinitely 
large. Since the interaction between the cracks can be neglected, the problem can be 
transformed into a single crack in an elastic half-plane. The numerical method 
described in the previous section can certainly be modified to solve for the stress 
intensity factor for this configuration, but we decide to use a simpler and more explicit 
approach. According to Tada (1985) (page 8.3a), the stress intensity factor can be 
expressed as 


K= 


2/, r a F(y/a) 
-fm™ ^1 - ( yfa ) 2 


'f&dy 

d 


(3-35) 


where F(x) = 1.3 - 0.3x 5/4 . This formula has an error less than 0.5%. 


Let 


K 2 r F(y/a) 

J, Vl-( y/a f 


a 


By introduce dimensionless variables 

s = y/a, z-y/d and t - a! d 


(3-36) 


(3-37) 


and dimensionless function 
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®(0 = 

Jo Vl-S 2 

(3-38) 

thus 

K 2 = -a® 2 (t) = —t® 2 (t) 
n Jt 

(3-39) 


The term K 2 must satisfy the third condition of crack initiation in the form 

ra — 2 — 2 

j a Kdy = aK 

It can be transfer to the following form by using the relation (3-38) 

£ zO 2 (z) dz-t 2 <^> 2 (t) (3-40) 

We can define another dimensionless function 


'V(t) = [z<b\z)dz 


and rewrite (3-40) as 


x ¥(t)lt = tQ 2 (t) 


Also according to (3-3 1) and (3-39) 


4 d 


l 0 - K =— 
n 


so two following equations are obvious existed 


it>‘( 0 = ~z VW /' = ■ n 


4 d* 


4 d’ 


(3-41) 


(3-42) 


(3-43) 


To achieve adequate precision, Gauss-Chebyshev quadrature must be employed to 
calculate the function O. First we extend the integration interval to (-1, 0) as 


(3-44) 



We can determining the ratio a/d from equation (3-42). Then, through either of 
equations in (3-43), we can determine the value of d*. This value of d* is found to be 
approximately 2.281, and the corresponding crack length a*=T.411 (when the initial 
strain profile is taken as a parabolic function). This result is also used as another check 
on our numerical calculations. The difference between the values just calculated and 
the extreme values of our previous numerical results is less than 0.4%. 

The function /<D 2 (r) and ¥(/)// is plotted in Fig.3-3. Note that, away from the 
origin, there is only one point at which these two curves intersect. Such a point 
happens to be the maximum point of function 'F(/)/f, which is the dimensionless form of 
the total energy released due to the crack formation. This property can easily be 
verified by the definition (3-38) and (3-41) of these two functions. The ratio t = a/d = 
0.6186 is such that the total energy released is maximized among all the other ratios. 
As a result, the penetration depth is minimized. 

One is naturally led to the question: what will happen when the maximum 
tensile stress exceeds the tensile strength of the material while the penetration depth d* 
is still much less than the minimum value 2.281? Such a situation can happen, for 
instance, when the half plane represents a very hot object, and the surface of which is 
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suddenly brought into contact with a very cold medium. The surface stress quickly 
rises to Ea(T x - T 0 ) (where a is the coefficient of thermal expansion and - T 0 is the 
temperature difference), but the penetration depth is initially very small because there is 
not enough time for the conduction of heat into the material. 

According to our theory, when d* is smaller than the lower limit value, there 
will not be enough energy available to open a crack, although the stress level is already 
high enough to break the material. As a result, the maximum tensile stress can rise 
above the tensile strength while the material retains its integrity. Now we can replace 
the stress condition with 

<7 x (°,y) = rf t f(y ,d ) ( 3 - 46 ) 

Then, equations (3-40) and (3-44) will change to 

* ! = — '4> ! (') = T3fr. (S- 47 ) 

n 4d*y 

When the stress condition is changed, the energy balance laws are also changed 
accordingly. In fact, the altered system will be the same as the equations in (3-42) 
except that the dimensionless loading depth d* has to be replaced with y 2 d*. As a 
result, the critical ratio a/d is the same and the minimum penetration depth becomes 
2.28 \ly 2 , which is still finite, although smaller than the original minimum depth. 

3.6 Experimental Evidence 

The present theory appears to be compatible with the existing experiment 
evidence. In Geyer’s experiment (Geyer et.al., 1982), a uniformly heated (at about 200 
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°C) glass plate was put into contact with dry ice (at -78 °C). Seconds after contact, a 
few cracks suddenly shot up in a dynamic manner. Since the paper did not report what 
kind of glass is used in the test, we assume that it was soda-lime glass. For this kind of 
glass, the tensile strength is typically 70 MPa (Bansal, 1986). The typical standard 
deviation in the glass tensile strength is about 20%. The Young's modulus reported in 
the paper is E = 69 GPa. The thermal expansion coefficient is 8.5 x KT^C. The initial 
tensile stress caused by the temperature difference is calculated to be of 60MPa, 
which is larger than the tensile strength. If the fracture energy G f is taken as 3.6 N/m 

as reported in the paper, the effective length of material l 0 is about 5 x 1 0“ 2 mm . With 
such a small reference length, only the part of the solution for small d* and big b* is 
relevant to this experiment. 

It is observed that cracks do not form immediately after the hot glass is in 
contact. Rather, a few seconds are usually needed. This can also be explained by our 
theory: the penetration depth needs to reach a certain length (a process that takes time) 
before the cracks can form. The average spacing is about 1 to 2 cm, which is 
sufficiently large for b* to be regarded as infinity. In this experiment, the maximum 
stress is larger than the tensile strength (X is about 2). An over stressed plate can be 
highly unstable. For instance, any disturbance supplying energy (such as kinetic 
energy) can make the energy sufficient for crack formation. If this condition is met, the 
crack could form suddenly and probably in a dynamic manner. Another possible cause, 
for dynamic crack growth, probably more important, is the fact that the tensile strength 
decreases with an increase of temperature. Although the fracture energy also decreases 
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with temperature, the rate of decrease is not as significant as the tensile strength. In the 
case of Geyer's test, the cold side of the glass plate is much stronger, therefore requires 
a higher level of stress to break the glass. However, once the cracks are formed and 
enter the zone of higher temperature, where the glass is less tough, then there is a 
surplus of crack driving force (that is, energy release rate). Such a surplus of driving 
force would certainly cause cracks to grow in an uncontrolled dynamic way. 

3.7 Additional Comments and Brief Summary 

Which d* should be used depends on the loading method? Imagining that both 
the maximum value of the initial strain and the loading depth grow with time, one 
should take the loading depth at the moment when the maximum strain reaches the 
tensile strength. It is this loading depth that determines the initial crack spacing and the 
initial crack length. However, if the loading depth is still smaller than the minimum 
value, then the critical penetration depth will be the depth that first satisfies the relation 
d* = 2.2%\/y 2 . 

Nemat-Nasser et al. (1978), (1979) and gave an estimate of the lower limit of 
crack spacing by energy consideration. The idea is that the fraction of total strain 
energy before cracking should be used to generate new surfaces. From which , using 
our notation, the minimum crack spacing in the form of b*>l/(y 2 K 0 ). K 0 is an 
empirical coefficient depending on the initial strain loading profile, the ratio a/d and a 
parameter 9 representing the fraction of strain energy to be used in creating new crack 
surfaces. However, neither a/d nor 6 can be determined in their analysis, as is admitted 
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by Nemat-Nasser (1979). In Geyer's paper, this factor is taken as 0.1. Therefore, their 
estimate of lower limit is 5 times larger than here. Bazant et al. (1979) also gave a 
similar estimate, and the value of a/d is considered to be at least 1.5, which is also 
larger than our result a/d = 0.62 (see Fig. 3-3). An adequate experiment is needed to 
verify these theoretical predictions. Glass is probably not the best choice for the test, 
because its effective length is so small. Crack sizes in the range of 1 to 10 times of /„ 
are very hard to observe by the naked eyes. Brittle materials with / 0 in the range of 
millimeters or even centimeters are preferable. 

For small d*, the initial crack length is of the same order as the effective length 
/ 0 . Note that / 0 is of the same order as the process zone length, which is a zone that 
contains highly nonlinear deformation. When the crack length is of the same order as 
the process zone length, the nonlinear effects are no longer negligible. The crack 
initiation theory, which is mainly based on linear elastic fracture mechanics, should be 
only regarded as the first approximation. The problem of determine the initial crack 
spacing in a material with a large nonlinear process zone is to be studied in the chapter 
V. 

The crack initiation theory proposed herein is more general than either 
homogeneous isotropic materials or a half-plane problem. For example, for orthotropic 
or anisotropic material you may find a suitable representation for the components of 
stresses or the stress intensity factor (as a function of crack length, loading depth and 
crack spacing), then the problem can be solved in the same way as for isotropic 
material . For the purpose of applying the theory to more practical engineering 
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problems and also checking the correctness of the initiation theory, the crack initiation 
of concrete pavements is to be studied in the next chapter. 

The main contribution done in this chaper van be briefly summarized as follows: 

1. The initial spacing and initial stable equilibrium length of parallel equidistant 
cracks emanating from the surface of a brittle elastic half plane can be determined from 
three conditions: (1) The stress at the surface reaches a given strength limit. (2) After 
the initial cracks form, the energy release rate equals its given critical value. (3) The 
finite energy release due to the initial crack jump equals the energy needed to form the 
crack (according to the given fracture energy of the material or fracture toughness). 

2 The problem can be solved if the stress intensity factor as a function of 
loading depth, crack length and the crack spacing is known. The stress intensity factor 
can be solved using Cauchy's integral equation. For the limiting cases of infinite initial 
crack length and of infinite large crack spacing, the correct limiting values are 
approached, which provides a check for the accuracy of the numerical solutions. 

3. The results of analysis compares favorably with available experimental 
evidence on thermal cracks in glass. However, this comparison is valid only in the 
range of very large initial crack spacing. For a complete check of the validity of the 
solution, further experimental studies are needed. 

3.8 Outline of a Proposed Experiment 

For a complete check of the validity of the solution, further experimental 
studies are needed. The following experiment may be conducted in a laboratoiy. 
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The brittle material suggested for experiment is gypsum, which is cheap and its 
characteristic length lo is suitable for a limit sized spacemen simulating a half-plane. 
Glass’s lo is too small the naked eye may not observe crack forming. On contrary 
concrete’s lo is too large, required size of specimen would be beyond our laboratory 
ability. 

It is important to evaluate the material fracture properties. Beam specimens can 
be used for three-point bending tests. The dimensions of the beam specimen could be 
12 inches long, 3 inches high, and 1 inch thick. During three-point beam test, the 
applied load, the CMOD (crack mouth opening displacement), and the load-point 
deflection (5) should be monitored and stored in a digital form. The rate of loading is 
controlled by a constant increment of CMOD. The suggested loading rate is 0.01 
in./min. The typical load-load line deflection (P - 8) and load-CMOD (P - CMOD) 
curves can be obtained from the test data. From these two curves. Young’s modulus 
and fracture energy (Gf) can be evaluated. The Young’s modulus value calculated 
from the load-CMOD curves are usually higher than those obtained from the load- 
deflection curves, as reported by other researchers. Fracture energy is defined as the 
area under the load-line deflection curve divided by the initial uncracked ligament area. 
Although the area under the load-CMOD curve does not have a direct physical 
meaning, the fracture energy calculated from it is comparable with that obtained from 
load-deflection curves. Since CMOD value is less sensitive to support settlements, the 
P-CMOD curve should be used to determine the fracture energy. 
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For determining the tensile strength , splitting tension tests can be conducted, 
which is carried out on a standard cylinder, tested on its side in diametral compression. 
The dimensions of the cylinder specimen could be 3 inches in diameter and 6 inches 
long. The test could be performed under displacement control at a rate of 0.005 in./min 
or under load control at a rate of 50 to 150 lb/in. 2 /min. The tensile strength then can be 
IP 

calculated as / = , where P is the applied compressive load, L the cylinder 

nLD 

length, D the cylinder diameter. 

To observe the phenomena of crack initiation due to drying shrinkage the 
following experiment may be conducted in some laboratory. The top opened rectangle 
box, which inside dimension may be three feet width, one inch thick and two feet high 
(the final size should depend on the fracture material properties), can be made by 
concrete or some hard material because the deformation in width direction is not 
allowed. Two side panel in thickness direction should be greased, it would allow free 
movement along these panel surfaces. The moisture measure meters may be installed in 
the middle section distributed from top to bottom of box. The prepared gypsum paste 
then can fully fill the box. Keep the top of box open to evaporate water from this top 
surface. 

The moisture record every certain period of time until the shrinkage cracks 
appear. Measure the crack spacing and depth and draw a moisture distribution curve, 
it can be modified as a drying profile. The results from the experiment and from the 
initiation theory can be compared and reported. 
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(C) 


Fig. 3-1 (a) Geometry definition of parallel crack system; 

(b) Initial strain profile; 

(c) Unit cell of width b. 





CHAPTER IV 


CRACK INITIATION IN CONCRETE PAVEMENTS 

It has been known that joints must be placed in Portland cement concrete pavement 
to control cracking that is caused by environmental changes. By providing joints, cracks 
either do not develop or develop only in a controlled manner related to the location of the 
joint. The cost of construction increases and the ride quality decreases if the spacing 
between the joints are short. Selection of joint spacing is motivated by above requirement. 
Many states have developed standards for joint spacing through experience. No rational 
procedure exists to determine the “correct” joints spacing. However, due to paucity of 
general understanding of the crack spacing phenomena, the best guideline, at the present 
time, is still the local road service record or certain rule of thumb. For instance, it is 
advised that the joints spacing (in feet) should not greatly exceed the slab thickness 
(in inches). It is therefore important to understand the basic physical mechanism of 
pavement crack spacing. 

For asphalt concrete pavement, thermal cracking is one of the main causes of quality 
deterioration. Each seasonal cycle brings about new cracks in the surface course, and the 
extent of cracking is usually measured by the cracking index, which is the length of cracks in 
each unit area. The theory proposed in this paper is relevant when pavement surface is free 
of cracks, or when the cracks are still very sparse. 
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This dissertation aims to establish a simple mechanical model of concrete pavement 
that can be used to predict the spacing in pavement due to temperature differential. Crack 
spacing due to moisture change can be discussed in the same way if the coefficient of 
expansion due to moisture change is known for the material. The main purpose of this 
study is to identify the basic laws that govern the cracking phenomena in pavement, and 
using these basic laws to relate causes to their consequences. With these knowledge, it may 
improve our capability to extrapolate the experimental findings in one locality to another 
and thus reduce the cost and time required for repeated field tests on different geographical 
locations. The model proposed hopefully can serve as a theoretical foundation for 
pavement joint design. 

Cracking in pavement is an example of a special class of mechanical problems in 
which a system of cracks suddenly occurs on a smooth surface under the action of 
temperature and moisture change. The problem of crack initiation from the surface of an 
elastic half-plane has been studied in the previous chapter. However, the problem of crack 
initiation in pavement has its own rights because of the following reasons. First, the 
pavement is very different from a homogeneous half-plane, the material property is 
dramatically different between the pavement and its foundation (or subbase). Second, the 
pavement material is not perfectly brittle so linear elastic fracture mechanics cannot be 
applied. As a result, a new modeling must be devised to take advantages of the structures 
main feature, and the crack initiation theory must also be modified to take into account of 
the influence of large process zone. 
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4.1 Mechanical Modeling of Pavement 

The pavement is structurally modeled as an elastic plate sitting on a Winkler 
foundation(Fig.4-la). To focus on the main issue, only transverse cracking will be 
considered in this paper. Therefore, we can take a slice of plate of unit width together with 
its supporting foundation in the longitudinal direction from the plate. Since we are only 
interested in the average behavior of cracking spacing in pavement, the crack spacing 21 and 
crack depth a are assumed to be uniform. Following the treatment of Okamura (1973, 
1975) and also that ofRice and Levi (1972), the effect of cracks in a beam is represented by 
the increase in its compliance. Let A and 6 1 be the total elongation and rotation of an 
elastic beam with an edge crack at the center (Fig.4-2) loaded by bending moment M 
(moment per unit width) and axial tension N (force per unit width) at both ends. M and N 
are reaction forces from adjacent plates. The positive M is defined as it would open the 
crack ( top of beam in tension). The total deformation can be expressed as the sum of the 
deformation of the beam without a crack and the additional ones due to the crack, that is, 

A, = A 0 + A c and 0 t = 0 O + 0 C (4-1) 

where subscript "0" denotes deformations of the beam without a crack, and ‘fc” additional 
deformations due to a crack. These deformations also are corresponding to the whole beam, 
which are the sum of both ends deformation. The additional deformations can be related to 
these nominal forces as: 


0 C — Xmm M Xmn N , A c - Xnm Ai + Xm N 


(4-2) 
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where X represents the increments of the compliance functions caused by the existence of a 
crack, namely, Xmm — additional rotation due to an unit bending moment, Xm = additional 
elongation due to an unit axial tension force, Xun = Xiu (due to reciprocity of linear 
elasticity) = elongation or rotation due to a unit value of M or N respectively. These 
additional compliances would be defined in section 4. 

The basic equation of lateral equilibrium of plate resting on Winkler foundation can 
be written as 

D^r+kM = 0 (4-3) 

d x 

where D = Ebh 3 /\2( 1-v 2 ) = plate stiffness and K = the coefficient of subgrade vertical 
reaction; E =Young's modulus of pavement, v = Poisson's ratio of pavement, b = plate’s 
width, and h = plate’s thickness. This equation is the same as that for an elastic beam in the 
plane strain condition. The external load need not be considered because only thermal effect 
is of concern. Also because the self weight of pavement only produce uniform downward 
movement and increase the magnitude of the contact stress, thus, the flexure of the beam is 
unaffected by effects of self-weight. We need not to consider it either. 

The temperature vary through the thickness of pavement but not in planes 
parallel to the surface of the pavement. For simplicity the temperature distribution 
through the thickness of the pavement is assumed to be linear. This assumption is 
without significantly affecting the accuracy of the computations., as Armaghani, Larsen 
and Smith (1987) pointed. 
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The pavement is stressed by temperature differential, with T, = the temperature 
differential at the top of the beam and T b = the temperature differential at the bottom. 
Denote 


Ea Tt-Tbh 2 „ EaT t +T bl _ 
Mt ~ — r Nt = «- 


1- v 2 


1-v 


(4-4) 


as the moment and tension that the beam would experience if the corresponding 
deformation is restrained, a = thermal expansion coefficient of pavements. 


4.2 Structural Analysis 

The structure can be solved by the compatibility condition at the cracked section. 
Due to symmetry, only one half of the beam needs to be considered (Fig.4-lc). The elastic 
solution of equation (4-3), found from Selvadurai’s book (1979), pp.84, is 

w(x) = e /Uc ('C/C 0 s>^ + C 2 sinAx^ + e'' ,J 'rC 3 C 0 s^x+c^sinAr / ) (4-5) 

where 'k 4 ~K I4D, X x is a dimensionless parameter. The constants Q, (i=l,2,3,4) are to be 
determined by the boundary conditions of the beam. 

From the homogenous solution (4-5) the slope , bending moment and shear force of 
the beam are given by 

= M(x) = D*^-M r - K(x) = -D*^f (4-6) 

ax ax dx 

Boundary conditions 0(0) = 0, V(0) = 0 and F(I) = 0 yield 

C, = C 3 ; C 2 = — C 4 ; C, = QC 2 


Then the solution can be written in the forms: 

0 O =2 9(1) = 4 C 2 W(Xl), M(l) = 4DC2l 1 ®(M)~M T 


(4-7) 


Where O and 'F are dimensionless functions given as 
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The function Q is plotted in Fig.4-3. As can be seen, when the nondimensional length A1 
becomes larger about 1.8, Q starts to fluctuate and eventually approaches to 4. However, 
this behavior is not desirable from physical point of view, because it amounts to say that the 
bending compliance is bounded from above, while in reality the bending compliance must 
increase with beam length. 

The assumption that the reaction from the foundation is proportional to deflection is 
reasonable when the beam deflection is downward. It becomes meaningless if the deflection 
is upward. A small portion of upward deflection is acceptable because the self-weight 
always produce uniform downward movement that keeps the beam in contact with the 
foundation, but if there are significant uplifting of the beam, then the compliance formula 
must be modified. The simplest remedy is to assume that the stiffness of the foundation can 
be neglected, one can use the simple beam theory to obtain 
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(4-13) 


2 6{l) 21 

M(l) D 

this compliance formula is proportional to the beam length. It will be seen later that (4-12) 
and (4-13) yields very close solutions for a small crack spacing, (4-13) is more reasonably 
behaved if the crack spacing is large. Therefore, the compliance function defined in (4-13) is 
preferable. 

The elongation is related to the elastic normal force and thermal force through the 
following relation 

\-Cw(N + N r ), OM 4 )- 

Now the total deformations at the ends of beam , can be expressed as 

9 t ~ Cum (M + Mt) + (Amu M + Amn N ) (4-15) 

A, = Cm(N + Nt) + (AnmM+AnnN) (4-16) 

The final equation for solving the bending moment and axial force depends on 
additional assumptions regarding the constraint in the longitudinal direction. 

Case I. The full constraint at the ends of beam : 9 t = 0 and A ( = 0. 

This case reflects that for an infinite long pavement , the ends of any unit beam are the 
surface of symmetry, at which above constrains conditions established. Through Eqs.(4-15) 
and (4-16), the bending moment and axial tension force can be solved from a set of coupled 
equations: 

(4-17) 


(Cmm +Amm)M + Amn N +CmmM t -0 
AnmM + (Cm + Am)N + Cm N T = 0 


(4-18) 
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As we will see later, the crack spacing will be very sensitive to the thermal contraction, 
which is presented in the equation in terms of Nt. 

Case II. Constraint only on rotation but not in axial direction: Q t - 0. 

If there is no bonding or friction between the pavement and its subgrade, then the beam can 
expand freely in axial direction and the axial force N is zero consequently. We obtain an 
equation to determine the total moment as 

(Cmm + Cmm M t = 0 ( 4 - 19 ) 

Case HI. The same constraint condition as Case H and plus friction force. 

If there is friction between the pavement and its subgrade, then during contraction 
the friction would act against the relative movement to prevent the contraction. Across the 
center of the beam the frictional force must change sign because of symmetry. The axial 
force at the end of the beam, denoted as Nf must be the total fractional force on one half of 
the beam span. Equation (4-17) must be modified as 

(Cmm + — Cmm M t — Xmn N f ( 4 - 20 ) 

The simple formula of friction force (lb.) proposed by Friberg (1954) was 
f (l - x)F FI 2 

N f = J - — c tx. - for a slab with length 21 inches, concrete weight of 1/12 lb. per 

cu. in., and average frictional coefficient F. According to his paper, the range of the ratios 
between elastic modulus and friction coefficient might be from 10 6 /0.33, to 4.5xl0 6 /1.5. For 
the typical concrete system, with E = 4.2x1 0 6 psi, the friction coefficient F should be taken 
as 1.4, thus the friction force N { = 1?/120. 
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The frictional force acts to enhance the effect of thermal loading in the pavement. 
This is in agreement with the engineering experience that the joint spacing must decrease for 
an increase of frictional coefficient. . However, according to our calculation, the effect of 
friction force on the crack spacing is practically negligible. 

Subbases used many years ago were soft, loose and unbounded, such as sand, 
gravel and clay. Today stabilized agents, such as cement, lime and asphalt, are added to 
subbase materials, which lead to a strong bonding between the pavement and subbase in the 
tangential direction. As demonstrated by Wesevich (1987), the horizontal force is basically 
proportional to displacement for concrete slab placed on cement or lime treated clay 
subbases, suggesting that the interaction may be modeled as bonding rather than friction. 
For further increase of horizontal displacement, the resisting force approaches a constant. If 
that is the case, equation (4-18) should be used instead of the model to be proposed as 
follows. However, as will be seen later, neither bonding nor frictional forces seem to be 
important for determining the crack spacing in the pavement. 

Case 4. The pavement is bonded to the subbase and full constrain (as Case I) at the 

ends. 

If the pavement is assumed to be bonded to the subbase and the tangential force at 
the bottom is assumed to be proportional to the horizontal displacement, then the 
proportional constant, denoted as fa and representing another elastic modulus of the 
properties of the subbase, must also be given. The equilibrium equation for a slice of 
pavement (of unit width) in the axial direction can be expressed as 


The solution that satisfies the condition 7/(0) = 0 can be written as 

u(x) = Csinh/cc, N(x) = C— ^-cosh/cc (4-22) 

where // = (\-v 2 )fa IEh and C is an arbitrary constant. The compliance function at the end 
of the beam is obtained as 

(4 - a > 

It is straightforward to show that when fa becomes vanishingly small, the compliance 
function increases and approaches the limit value Cm ~2 l(\-v 2 )/Eh. Therefore, bonding 
compliance reduces the axial compliance only, as our intuition would expect. The total axial 
force and bending moment must be solved by Eqs. (4-17) and (4-18) with Cm being 
calculated according to (4-23). The moment contribution of this bonding force is neglected 
for the benefits of simplicity of the formulation. 

Before we proceed further, it is important to note that in the above discussion, the 
thermal load is calculated on the assumption that the subbase dose not contract or expand 
with temperature. If the thermal deformation of subbase is also considered, then the thermal 
loads defined by (4-4) must be modified accordingly. 

4.3 Stress Intensity Factors and Additional Compliance Functions 

The stress intensity factors at crack tip caused by an axial force and bending 
moment, according to Tada’s handbook (1985), can be represented respectively as: 


K * = 7h kJe) ’ Ku= ^h ku(e) 
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( 4 - 24 ) 


where e = alh = relative crack length and 


fc w(*) = ^2tan 


— 0.752 + 2. 02x + 0.37(1- sin —V 

2 _ 


cos- 


;rr 


( 4 - 25 ) 


— 0.923 + 0.199(1 -sin—) 4 
2 


=4 


M*) = 6 -t 2t an 


cos- 




( 4 - 26 ) 


The total stress intensity factor AT = Ajv + A*/. 

The elastic strain energy of structure due to cracking (also substituting the eqs (4-2) 
into expression) is 




( 4 - 27 ) 


The energy release rate in terms of stress intensity factors by Irwin’s formula can be 


expressed as 

G = (K„ + K u f(\-S)IE 

Also g = - = -N 2 ^^ + NM^^ + -M 2 ^^ 

da 2 da da 2 da 


( 4 - 28 ) 


( 4 - 29 ) 


By comparing above two expressions of G the following relations established between 
compliance functions and stress intensity factors as: 
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dk m 2(1 - v^) , 2 
Ja Eh N 


<&NM _ 

da 


2(1- v^) 
Eh 2 
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^ ML " 2(1 - 3 l/) ^ ( 4 - 30 ) 
c/a £// 3 w V ' 


By integrating them the compliance functions can be expressed as 

(4-31) 

r\ h. 


where i, j = M or N, tj = 0 for Am, 77 = 1 for /W and tj = 2 for The compliance 
function can be calculated according to these formulas using numerical integration method. 
In fact, the expression for and Am can be obtained directly from the fracture 
mechanics handbook compiled by Tada et al. (1985), only Am needs to be calibrated 
according to (4-3 1) using numerical integration, see Appendix 2. 


4.4 Crack Initiation Theory 

The previous structural analysis allows us to calculate the internal bending moment 
and axial force once the crack spacing 2/ and crack depth a are known. If there are no 
initial cracks in pavement, the internal moment and axial force can also be obtained using 
the method discussed in the previous section by letting the compliance of the cracked 
section being zero. As the thermal loading continues to increase, the tensile strength of the 
pavement is reached, and the necessary condition for crack initiation is satisfied. However, 
the strength criterion cannot determine what is the average crack spacing, nor can it 
determine the average depth of the cracks. 

The energy release rate can be related to the net (total) stress intensity factor K by 
the Irwin formula as K 2 ( 1 - v 2 ) / E , the total stress intensity factor can be expressed as 


l 
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( 4 - 32 ) 


K = K„+K u = -^k„(e) + -^Lk u (e) 

where N and M is determined by the method discussed in the previous section. 
Although the function kn and k M increases with relative crack length aJh, N and M 
decreases with aJh because if the section is totally severed N and M becomes zero due to 
loss of constraint. Once cracks start to gain enough lengths, the crack will stop growing 
when Griffith condition G = Gf is satisfied. During the crack growth, the total energy 
needed to form new crack surfaces, which is aGf for each crack, must be equal to the 
energy released by the structure, which is the shaded area under the curve in Fig. 4-3. This 
is the very condition that determines the crack spacing. 

The process of crack initiation consists a transition from the pre-initiation state, at 
which the strength criterion is satisfied, to the post-initiation state, at which the Griffith 
criterion is satisfied. More specifically, the proposed initiation theory in Chapter HI can be 
restated as 

(1) The maximum tensile stress must be equal to the tensile strength of the 
pavement at the pre-initiation state; 

(2) The energy release rate of the structure must equal to the pavement toughness 
G c at the post-initiation state; 

(3) The total energy must conserve during the initiation process. 

These three conditions play different roles in the solution procedure. The first 
condition determines the critical value of temperature differential, but it is not related to the 
cracked configuration. The temperature profile is considered fixed during the initiation 
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process. The rest of the two conditions must be solved together to determine the crack 
spacing and crack depth for the post-initiation state. 

Mathematically, with the temperature profile determined by the first condition, the 
axial force jV and the moment Mean be solved as functions of half crack spacing / and crack 
depth a. I and a can be determined by the second and the third condition, which are 

G(a) = Qf, Yl(0)-U(a) = aGf (4-33) 

where P is the elastic strain energy of the structure. According to the definition of the 
energy release rate, the third condition can also be expressed as 

\ a 0 G(a')da' = aG c ( 4 - 34 ) 

The integral usually must be evaluated by numerical quadrature. 

In the previous discussion, the material is assumed to be perfectly brittle. In terms 
of fracture mechanics analysis, brittleness means that the nonlinear process zone around the 
crack tip is negligible compared to characteristic structural dimension. It is widely known 
that Portland cement concrete is not a perfectly brittle material, its nonlinear process zone is 
often of several inches in size (Bazant, 1986a). To take the nonlinear behavior of the 
material into consideration, one must perform fracture mechanics analysis according to, for 
instance, the cohesive crack model (Hillerborg, 1976) with the softening stress described by 
the crack opening displacement, or equivalently, a crack band model (Bazant and Oh, 
1983). The crack initiation theory based on nonlinear fracture mechanics theory, although 
extremely interesting and will be pursued in a future study, is too complex for the purpose 
of this study. 
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To obtain a simpler model, one should be cognizant of the difference between the 
fracture energies obtained using different methods of calibration. Based on the cohesive 
crack model, one can obtain fracture energy based on the work of fracture, as is shown 
schematically in Fig.4-4. The total work FFioui done by the load P during a displacement 
controlled experiment is assumed to be totally transformed into creating new crack surfaces, 
which is equal to 2B(h-a), where B = the beam width, therefore the total fracture energy G/ 
= WtotJ2B(h-a). For Portland cement concrete, it is typically 80 ~ 120 N/m. 

There are many other methods to calibrate the toughness of concrete on the basis of 
maximum load capacity, among them we mention the size-effect approach proposed by 
Bazant and Pfeiffer (1987) and the two-parameter model proposed by Jenq and Shah 
(1985). These methods are based on the peak load value of the experimental measurement. 
The fracture energy so obtained would be denoted as Gf, which is usually in the range of 
30-60 N/m. 

The difference in the values of fracture energy can be understood from the point of 
view of the cohesive crack model. The maximum load that a structure can sustain, 
according to the cohesive crack model, depends on the slope of the softening curve, not the 
total area under the softening curve. This is because that the crack tip opening displacement 
at the peak load is less than its threshold value w c . For concrete slab with thickness around 
10 inches, the crack tip opening displacement at peak load is typically less wJ2. In other 
words, the shape of the softening curve that corresponds to values larger than wJ3 is 
irrelevant for the purpose of maximum load determination. In other words, the softening 
curve for the purpose of maximum load can be effectively replaced with a linear softening 
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law with its slope corresponds to certain average slope of the original nonlinear softening 
law in its small crack opening displacement range. The total area under this linear softening 
law is denoted as Gf . As shown by experimental studies, Gf is typically 1/2 ~ 1/3 of G/, 
which is the total area of the original nonlinear softening law, as is shown schematically in 
Fig.4-4. 

The second condition of crack initiation is concerned with the load bearing capacity 
of the structure in the post-initiation state, therefore it seems appropriate to use Gf as a 
measure of material toughness. The third condition is a statement of energy conservancy, 
therefore Gf appears to be more relevant as a measure of material toughness. It should be 
noted that the difference between the two measure is size-dependent, and Gf would 

approach Gf if the characteristic size of structure becomes very large. Thus, the crack 
initiation theory given previously is correct if the material is perfectly brittle. But for 
pavement, the difference between Gf and G/ must be reflected in the crack initiation 

theory. To this end, the second condition of the crack initiation theory may be written, 
mathematically, as 

G = Gf = PGf, P = Gf /G f (4-35) 

where the typical value of (3 ranges from 1/2 to 1/3. The third condition remains the same. 
In this way, the nonlinearity of the pavement material can be approximately accounted. As 
will be seen later, the modification dose not have significant influence on the final crack 
spacing and crack depth. The nonlinear effect of the material, therefore, does not seem to 
be consequential by this approximation. 
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4.5 Numerical Method 

It is useful to use a non-dimensional form to understanding the dependence of the 
solution on various material properties. For this purpose, the following non-dimensional 
nominal stresses can be introduced 


= 


jV_ 

f h' 


M 


6M 

f.h 2 


(4-36) 


where f is the tensile strength of the material. With these notations the first condition of 
crack initiation theory can be simply expressed as 

+ = J (4-37) 

while the second condition becomes 

+ ^ (4-38) 

where lo = E Gf )f? = the effective length of the pavement. The third condition 
becomes 

+ de’ = | (4-39) 

It should be emphasized that, although not explicitly written in (4-38) and (4-39), ctjv and 
a m depends on the relative crack depth as well as crack spacing. 

Also the equations for solving bending moment and axial force can be written in 
nondimensional form. By defining the nominal thermal bending stress as = 6M T / h 2 f t 
and nominal thermal normal stress as cf s = N r I h, the basic coupled equations (4-17) and 
(4-18) ( in Case I ) in the nondimensional form are 
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( Cmm + Xmm)&m ~ X\ in + Cmm — 0 (4-40) 

Xmn <j m!^ > + (Cm + Xm)&N + Cm °w = 0 (4-41) 

where the non-dimensional compliances due to a crack are defined as 

= \& k “ (e ' )de ' ' ( 4 - 42 ) 

Eh 

Xmn ~ Xmn~ — 2 \ 0 kn( e> ) kM( e ' )de' (4-43) 

Fh 

Xm ~ Xm~^ = 2j g kj/( e' )de' (4-44) 


and their empirical expression are presented in appendix 2. 

The non-dimensional compliance functions are defined as 
E 


Cnn ~ Cm 


l-f 


2X1 . h 2 E 2 Q(Xl) 

i *■ Ci " ~ c *"7rv - it - (4-4S) - 


Since the compliance function due to cracks is zero if the crack length is zero, the 
nominal stresses at pre-initiation state can be simply solved as cr M = - oj and 
cr N = -<yj . Therefore the first condition, the necessary condition for crack initiation, 
becomes 

° 7 v/ + ^=-l (4-46) 

which becomes a condition to determine the critical temperature deferential. When the 
distribution of the thermal stress is far from being linear, then (4-46) should be disregarded 
in favor of the condition expressed in (4-37). 

Combining (4-38) and (4-39), we obtain the following equation: 
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f|( (4-47) 

For a given reasonable value of crack spacing /, one can always determine e from (4-47). 
To prove this, one needs to know that the net stress intensity factor K is zero when e 
approaches 0 or 1. The first zero is due to the fact that stress intensity factor approaches 
zero as e 172 when e— >0. The second zero is due to the fact that the total stress intensity 
factor K also approaches zero when e approaches 1. This is because that the compliance 
functions due to crack approaches infinity as 1/(1- e) 2 , forcing the nominal stresses due to 
tension or bending approaches zero as (1- e) 2 Although the stress intensity factors due to 
unit load approaches infinity as 1/(1- ef a , thus the total stress intensity factor K must 
approach zero like (1- c) 172 when e approaches 1. Translated into the language of energy 
release rate, we know that G approaches zero like 1-e when e approaches 1. 

The actual structure of K is more complex. Depending on the tensile component of 
thermal stress On T , K can be negative for a significantly large crack length, with the negative 
portion of K is deliberately cut off see Fig.4-5. This phenomena is particularly obvious for 
large compressive axial thermal stress cjn T (negative values in the plot). Although K 
becomes positive again for even larger crack length, this region cannot be reached unless 
additional energy is available to assist the crack to jump through the valley of negative K. 
Although the energy release rate is still positive for negative K, this value does not represent 
the energy that can be released from the structure by crack propagation. For this reason, it 
is necessary to ensure, by means of numerical algorithm, that the solution of crack depth 
does not cross the valley of negative K values. 
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Based on these observations, we know that G must be an increasing function 
initially and becomes a decreasing function eventually. When G is an increasing function, its 
average value from 0 to e must be smaller than the current value. Since (3 is less than 1, the 
left hand side of (4-47) must be less than the right hand side of (4-47) for small e value. For 
e closes to 1, the average value of energy release rate must be a positive number while the 
energy release rate approaches zero, thus the left hand side of (4-47) must be larger than the 
right hand side of (4-47). Consequently, there must be a relative crack length e between 0 
and the boundary of negative-.^ value that satisfies (4-47). This property is essential to 
design a successful numerical algorithm to solve the unknown crack spacing and crack 
depth. Once the crack depth is solved as a function crack spacing, one can use the second 
condition (4-38) to determine the crack spacing. This is the main feature of the numerical 
method used in this study. 

4.6. General Behavior of the Model 

In the following numerical examples, the thickness of the concrete slab is taken as 
10 inches, which is a typical value in the pavement design. The material properties used in 
these calculation is E = 4.2* 1 0 s psi, v = 0. 1 8 and K = 400 lb/in 3 . With these parameter, X = 
0.0293 (1/in.). The parameter (3 is chosen to be 1/3. For typical Portland cement concrete 
the material length is about 7 inches, in which the toughness of the material must be defined 
by the work-of-fracture method. The tensile strength f and thermal expansion coefficient a 
need not to be specified at this point because these parameters do not enter into equations 
(4-38) or (4-39). They will determine the critical value of temperature differential. 


74 

The simplest case is the pure bending without axial constraint (case 13). The crack 
spacing and crack depth in inches are plotted in Fig. 4-6 as a function of material length lo, 
which is also measured in inches. The effect of different compliance functions is clearly 
demonstrated. If Cmm is defined by (4-12), that is, if the pavement is modeled as a beam 
resting on a Winkler foundation, then there is a upper limit for l 0 beyond which there is no 
solution for the crack initiation problem. This upper limit is obviously related to the upper 
bound of the compliance function, as explained in the previous discussion. On the other 
hand, the crack spacing is basically a linear function if Cmm is defined by (4-13), which 
appears more reasonable from physical point of view. When l 0 is small, both compliance 
functions yield very similar results. 

When the axial deformation is constrained (case I), crack opening is restrained by 
the compressive force it generates, the crack depth is smaller whereas crack spacing is 
larger for the same lo ., as can be seen from Fig.4-7. The crack spacing approaches infinity 
when k. approaches certain finite value. One may argue that the extreme sensitivity of the 
crack spacing when lo is close to its upper limit value is not reasonable from the physical 
point of the view. For small lo. values, these two compliance functions again yield very 
similar results. Therefore, only the compliance function defined in (4-13) will be used in 
further discussion. 

The bonding between the pavement and foundation is plotted in Fig.4-8 with the l 0 . 
fixed 3 and 7 in., respectively. Adequate data for the k h value are difficult to find in the open 
literature. If we assume kh = K, then ji w 3xl0' 3 (in.* 1 ). At this magnitude, the tangential 
bonding is not particularly important if l 0 is 7 inches. It can be important for larger value of 
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Iq though. For smaller l 0 value this bonding effect can be neglected for all practical 
purposes. 

The effect of (3 is plotted in Fig.4-9 with no tangential bonding but with axial 
deformation constrained. Although it may come as a surprise, the influence of 3 is not 
important at all. The most substantial influence on crack spacing, in addition to the material 
length Jo, comes from the thermal stress distribution. As shown in Fig. 4-10, if there is a 
substantial tension component in the thermal loading, then the spacing is significantly 
reduced. If there is a compressive component, then the spacing is dramatically increased. 

It is instructive to know what will happen when there is only uniform temperature 
drop when the axial deformation is constrained. Since there is no bending, the beam would 
simply break into segments. The average length of these segments would be the crack 
spacing in this limiting case. Because the beam is totally severed, the crack depth a = h, the 
beam depth. Therefore the second condition is nullified, and only the first and the third 
conditions remain. The third condition can be written as 2 IJv^HE' = hGf The conversion 
must occur at a =/,. As a result, one obtains that the minimum crack spacing 21 = 2 1 0 . It 
can be seen from this analysis that lo is simply the ratio of specific elastic energy to fracture 
energy, which is why lo is one of the main factors that determine the crack spacing. This 
minimum value of crack spacing was also obtained in the previous study of spacing of 
parallel cracks initiated from the surface of an elastic half plane. 

For more general loading configuration, the elastic energy that can be stored, as 
well as the fracture energy needed, depend on the crack depth. This is why the second 
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condition of the initiation theory must be employed. However, the main picture is still the 
global energy conservancy. The more elastic energy is stored per unit length, the smaller 
the crack spacing will be; and vice versa. For instance, the elastic energy per unit length 
would be smaller if the restraint by the foundation is considered. That is why the spacing 
predicted with a Winkler foundation in effect is generally larger than the spacing without a 
foundation. For the same reason, the axial constraint leads to larger crack spacing because 
the deformation is constrained, thus less elastic energy is stored in the beam per unit length. 

The effect of axial thermal loading can also be explained. As can be seen from Fig.4-5, 
with a uniform contraction in thermal loading, the overall K value(and thus the total energy 
that can be released from the structure) is increased, thus smaller crack spacing. On the 
other hand, with uniform thermal expansion, the energy that can be released from the 
structure becomes smaller, thus larger crack spacing becomes necessary. 

4.7 The Effect of Nonlinear Temperature Distribution 

As was shown by Armaghani et al. (1987) that the nonlinear temperature 
distribution caused by daily temperature fluctuation is not very significant. In other words, 
the thermal load calculated according to (4-4) is a good approximation if we are only 
concerned with the thermal stress caused by daily temperature changes. However, any 
pavement structures with adequate design should not develop cracks under normal daily 
temperature changes. It is very likely that the critical condition is met when unusual 
weather conditions occur. For instance, a rainstorm in a hot summer afternoon will bring 
drastically cool off the temperature in the top layer of the pavement, as shown in Fig.4-1 la. 
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The temperature distribution becomes very different from a linear distribution, thus the 
thermal loads calculated according to (4-4) is no longer appropriate. Asphalt concrete 
pavement invariably develops thermal cracks in northern American states during winter 
season when the temperature drops rapidly in the night. The temperature profile is 
schematically plotted in Fig.4- 1 lb. In these cases (4-4) must be modified. 

Let 7(z) be the temperature increment (that is, the measured temperature minus the 
reference temperature at which there is no thermal stress) along the pavement depth. The 
thermal stress is £a77(l-v). The critical temperature distribution is determined by the 
condition that the stress at the top surface is equal to the tensile strength of the material. 
The thermal bending and thermal tension can be calculated as 

K = T(z)zdl, N r = -^2- J* 1 \z)dz (448) 

1 - v J ~7i \- v J -%. 

The nondimensional nominal thermal stresses Gm T and a N T are, again, defined according to 
(4-36). The crack initiation condition on the nominal thermal stresses (4-46), which is 
established on the assumption of linear temperature profile, must be discarded. 

It is assumed that the elastic bending and tensile stresses can still be determined by 
(4-40) and (4-41) during crack propagation. This is because once cracks occur, only the 
bending and tension is the dominant deformation modes away from the cracked cross 
section, which are captured by our mechanical model. As a result, cracks are initiated with 
lower level of elastic energy per length stored in the pavement, which in turn causes larger 
crack spacing. This effect can be important when theoretical predictions are to be 
compared with field observations. 
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As a demonstrative example, let's estimate the crack spacing in the asphalt concrete 
pavement. In the state of Illinois, the crack spacing developed in the pavement, after first 
winter, is about 8 to 12 meters on the average. Similar value of the first year crack spacing 
can be also deducted from the data given by McLeod (in the discussion of the paper by 
Ruth et al. 1982). The asphalt concrete is a very complex material, with material properties 
vary strongly with temperature. Asphalt concrete becomes quite brittle under low 
temperature, especially under temperature below freezing point. Although the experimental 
data of material length 1 0 is not available, we take lo = 10 in. for asphalt concrete, which is 
believed to be a reasonable value, judging from the fact that Portland cement concrete 
typically has l 0 =7 in. Assume also that the critical thermal stress ct m t = 0.7 and g n t = 0 
during a rapid temperature drop in a winter night. Furthermore, it is assumed that the 
pavement thickness is 10 in. and the effect of foundation can be neglected (i.e. the elastic 
bending compliance is calculated according to (4-13) and tangential bonding between the 
pavement and its subgrade is considered). Under these assumptions, the calculated crack 
spacing is 21 = 388 in., which is approximately 10 meters and is very realistic in comparison 
with the field observation. If the nonlinear critical temperature yields Gm T - 0.5, then the 
crack spacing would become more than 800 in., or 20 meters, which is still reasonable. The 
actually values of these parameters must be determined by experiment. 

4.8 Discussion and Conclusions 

In this chapter, a simple mechanical model for predicting the crack spacing in 
concrete pavement is proposed. The pavement is modeled as a simple beam with or 
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without a foundation, and the influence of cracks in the pavement is simulated by the 
increase of the compliance. As a result, a simple yet effective analytic model is obtained. It 
is emphasized that the occurrence of cracking in pavement is determined by the strength 
criterion, and no attempt is made to improve the strength criterion. Our model is aimed at 
predict the average crack spacing and crack depth after the strength criterion is satisfied. 
Therefore, the proposed model complements the existing knowledge on the phenomena of 
pavement cracking. 

The main contribution can be summarized as follows: 

1. The crack spacing can be influenced by many different factors. The most 
important factors are the material length l 0 and the distribution of thermal stress. The effect 
of the foundation can be important only under extraordinary condition, such as perfect 
strong bonding. 

2. For nonlinear material, a distinction must be made among the concepts of 
fracture energy. The effective fracture energy, which is defined on the basis of load 
capacity, can be as little as 1/3 of the total fracture energy, which is defined by the work of 
fracture method. The effective fracture energy must be used with the second condition of 
crack initiation theory, while the total fracture energy must be used with the third condition. 

3. It is demonstrated that the most main mechanism that controls the crack spacing 
in the pavement is energy conservancy. If the elastic energy per length that can be released 
by crack propagation is large, then the crack spacing is small, with 2l 0 being the smallest 
possible value. This relation can probably be developed into an experimental procedure to 
determine the material length /<? directly 
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4. Thermal cracking usually occurs during extreme weather conditions, in which 
the thermal stress distribution can be far from linear. When such a thermal stress becomes 
critical, meaning that the maximum tensile stress reaches the tensile strength of the material, 
the crack will initiation under rather small stress level, causing larger crack spacing. 
































CHAPTER V 


INITIATION THEORY — FOR QUASI-BRITTLE MATERIALS 

In the previous discussion, the material is assumed to be perfectly brittle, so that 
the linear elastic fracture mechanics can be utilized to describe the fracture behavior. 
Although the nonlinear process zone is taken into account for the concrete material 
used in highway pavement, the treatment is only approximate. In reality no material 
can be regarded as perfectly brittle. Perfectly brittle material is an idealization of the 
so-called quasi-brittle material when the nonlinear process zone is very small compared 
to crack length. In quasi-brittle materials fracture is preceded by a fracture process 
zone. It mainly consists of microcracking or other material defects. Experiment 
observations of concrete indicate that within this zone, fracture is manifested in the 
form of aggregate debonding and overlapping micro-cracks, mainly parallel to each 
other, with some intact grain bridges between them (Mier, Rots and Bakker, 1991). For 
illustration purposes this region is portrayed in Fig 5- 1(a). The cohesive crack model 
describes the region as a process zone (or a fictitious crack) where the material can still 
transfer stresses which decreases with increasing crack opening. In the crack initiation 
problem, the initial crack length is zero, therefore the process zone length is always 
important to be included in the analysis. The theory of cohesive crack model to be 
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employed in this dissertation is a nonlinear fracture mechnics theory, which was 
proposed by Hillerborg et al. (1976) for concrete material. 

When a structure contains more than one crack, the question of how cracks 
interact must be answered. These cracks may grow simultaneously, or some cracks may 
close during certain stages of crack growth. In the framework of linear fracture 
mechanics, Bazant et al. (1977, 1978), Nemat-Nasser et al.(1978) studied the 
interaction of equidistant parallel surface cracks. The problem of cohesive crack model 
with multiple cracks does not seem to be studied before. 

In this chapter the basic concept of the cohesive crack model and basic 
conditions governing the cohesive crack initiation will be introduced first. The crack 
initiation is more complex in cohesive crack model because it may involve the crack 
bifurcation, crack growth and crack closure. In particular, we propose the concept of 
unloading of a cohesive crack model to describe the crack closure. The half-plane with 
a system of alternative surface cracks is still formulated as a singular integral equation. 
The basic behavior of cohesive crack model with alternative crack lengths are 
discussed. To study the solution behavior, the rate equations (rate form of the basic 
equations) are developed. Based on the rate equations, the critical condition of CCM 
with interactive cracks can be reduced to an eigenvalue problem. Through it the 
maximum load and the bifurcation point can be directly found. The post-critical 
behavior is also discussed. It is very interesting that the lower limit of crack spacing in 
cohesive crack model is same as the one obtained for brittle materials. 



5.1. Basic Conditions and Definitions 
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In the theory of the cohesive crack model, it is assumed that fracture under 
monotonically increasing mode-I loading occurs when the maximum (tensile) principal 
stress reaches the tensile strength of the material f t . It is further assumed that fracture 
is localized in the so-called process zone such that there is no energy dissipation in the 
bulk of the structure. The process zone (Fig.5-l(b)) is modeled by a displacement 
discontinuity (a crack opening displacement) with the proviso that the faces of the 
discontinuity are capable of transmitting certain cohesive stresses, less than f t , such 
that a- cr(w) with<r(0) = /,, and cr(vt') > 0, where cr(w) describes the tensile 
softening behavior. In practice, cr(w) was approximated by a linear or bi-linear relation 
and was assumed to vanish when the crack opening displacement w reached a certain 
critical value w c . The fracture response in this model is characterized by fracture energy 
G f defined as the energy required to open a unit crack fully, or in other words as the 
area under the strain-softening curve between w = 0 and w = w c . In particular, we 
employ the linear softening law in which the cohesive stress in the crack surface can be 
expresses as 

°-= /,(!-— ), w<w c (5-1) 

where w c = 2 G f / f is the threshold value of crack opening beyond which a= 0. 

The stress distribution on the crack surface should include this cohesive stress 
that opposites to thermal loading stress and push the crack closure. Now the total 
stress distribution on the crack surface can be expressed as 


-F(0 + gO*(0) 
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(5-2) 


where F(£) = y/</(£/*0 is the thermal loading stress due to temperature 
penetration. For satisfying the strength condition the loading parameter y > 1 is added 
here because this thermal loading stress must overcome the opposite effect of the 
cohesive stress. Actually the thermal loading have two parameters: the maximum 
loading stress y and the loading depth d. They may grow simultaneously with time or 
only one of them grows. For simplicity, we choose y as loading parameter because F is 
the explicit function of y. 

In the previous work, we proposed the three governing conditions for brittle 
elastic solid. The first one, the strength condition, can still be used in CCM for the 
initiation of process zone. The crack initiation in quasi-brittle structure is a continue 
process in which the stress equilibrium condition will be the governing condition 
instead of the strength condition. It should be pointed out that the process zone 
initiation is not the same as the crack initiation. The second one is replaced by the 
condition of finite stress, which is discussed below. The third one is no longer 
applicable for the cohesive crack model. 

The stress in the cohesive crack model is finite everywhere, thus the total stress 
intensity factor at crack tip, which is the sum of stress intensity factor due to the 
external load and due to the cohesive stress, must be equal to zero: 

K(a) = 0 (5-3) 

This condition is initially introduced by Barenblatt (1959, 1962), and is the fundamental 
feature of any cohesive crack models. Compared to our previous study with linear 
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fracture mechanics, equation (5-3) corresponds to the Griffith theorem, that is, the 
energy release rate due to incremental crack growth must be equal to the fracture 
energy of the material. Furthermore, equation (5-3) implies that the total energy 
release rate is zero, because the energy release rate is proportional to the square of the 
total stress intensity factor. The condition that the energy release rate is equal to zero 
means that during crack growth (or to be more precise, process zone growth) all 
energy released from the deformation energy in the bulk of the structure is transformed, 
completely, to the formation of crack surfaces or overcome the work of the cohesive 
stress. (Although in linear elastic fracture mechanics the energy is also require to be 
converted into surface energy, the process of conversion is not included in the analysis. 
In cohesive crack model, the conversion of energy is realized by the work of the 
cohesive stress, which is a function of crack opening displacement, in the process zone 
ahead of the stress-free crack tip.) 

In our previous analysis, the Griffith theorem is not satisfied during the crack 
jump, therefore this condition must be supplemented by the condition that the total 
energy released during the crack jump must totally converted into the surface energy of 
newly created cracks, which is an independent condition. In the case of cohesive crack 
model, however, (5-3) is always satisfied. It implies that the third condition is no 
longer an independent condition. This difference in the energy consideration makes the 
initiation theory based on cohesive crack model fundamentally different from the 
initiation theory based on linear elastic fracture mechanics. In the cohesive crack 
model, we need to extend our understanding and concept of crack initiation, and find a 
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new approach to describe the process of crack initiation. This is the main goal of this 
proposed research. 

For quasi-brittle materials, the cracking of the material is a gradual process, the 
crack opening involves a large process zone in which the microcracks exist. If there is 
still cohesive stress acting across the crack surfaces, the crack should not be regarded 
as completely opened. In other words, we define that a crack surface is fully created if 
the corresponding crack opening displacement exceeds the threshold value w c . Also we 
need to define the crack initiation of the cohesive crack model at the condition of 
crack opening. In other words, if crack surface be fully occupied by cohesive zone, the 
crack have not initiated yet. 

5.2 Cohesive Crack Model with Interactive Cracks 

So far we have only considered the case where all cracks grow in equal rate. If 
the load can only increase to certain level which has not sufficient energy for all cracks 
to initiate, we want to know whether all these cracks stop growing or some cracks stop 
or even close in favor of other cracks to grow in a faster rate, resulting in larger crack 
spacing. To study this type of crack behavior, it is necessary to study a cohesive crack 
model with interactive cracks. In other word, we must study the cohesive crack model 
in which different cracks can have different lengths. These cracks may grow 
simultaneously, or some cracks may close in favor of other's growth. For simplicity we 
assume the every another crack has the same length, see Fig 5-2(a). 
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When a structure contains more than one crack, the question of how cracks 
interact must be answered. In the framework of linear fracture mechanics, Bazant 
(1977, 1979), Nemat-Nasser et al. (1978) have studied the alternative crack growth by 
differentiating the stress intensity factors with respect to crack lengths. However, the 
cohesive crack model with multiple interactive cracks has not been studied. The 
method mentioned above by other authors is no longer useful in cohesive crack model 
because the stress intensity factors equal zero. Two basic conditions: the stress 
equilibrium condition and crack tip conditions are the basic tools to probe the 
phenomena and help us to find out new condition in our problem. 

The crack initiation of cohesive crack model is a complicate problem that may 
involve crack growth, crack bifurcation and crack jump. 

Also we need define the crack unloading condition to discuss the phenomena of 
the crack closure. A cohesive crack is called in loading condition if its crack opening 
displacement increases, otherwise it is called in unloading condition. For simplicity it is 
assumed that all cracks are in loading condition in the preceding discussions. It is 
usually sufficient for CCM to study the load-deflection curve with a single growing 
cohesive crack. This is perhaps the main reason why there is little study about the 
unloading behavior of CCM. When there are more than one cracks, it is possible that 
some cracks are in loading condition while the other in unloading condition, even when 
the applied load increases. 

When a cohesive crack is in unloading condition, the cohesive stress does not 
follow function cr to increase when the crack opening displacement w decreases. The 
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material in the process zone is partially damaged, and the damage can not be reversed 
by crack closure. While the actual behavior of the unloading stress-displacement 
relation is a subject of experimental study, a simple relation is proposed as shown in 
Fig.5-l(b). Upon closing, the stress reduces linearly back to the origin. If the crack is 
reopened, the stress increases along the unloading line until the softening curve / is 
reached, then it decreases again following the softening curve a. This relation implies 
that the material can close perfectly to its original position. In reality the fractured 
surfaces are rough, the crack opening displacement may be unable to return to zero. 
One may need to use more sophisticated unloading relation to describe the immature 
crack closure. However, the precise nature of the unloading is beyond the my research. 

The cohesive stress at the crack unloading condition (in the linear assumption) 
is calculated as 

^>„)= ^ (5-4) 

w 

where w* denotes the crack displacement at the bifurcation point and w u of unloading 
cracks. This formulation insure the cohesive stresses linearly return to zero when the 
crack is finally closed. 

A loading crack implies a propagating crack. When a crack propagates, one 
requires that the total stress intensity factor be zero. The position of the process zone 
tip is unknown for a propagating crack, so the condition of zero stress intensity factor 
can be viewed as a condition to determine the crack tip position. On the other hand, if 
a crack is unloading, it no longer propagates. The process zone tip position of an 
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unloading crack is thus not a variable, and its correspondent crack tip equation is 

» 

redundant. Consequently, whenever a crack is in unloading condition, the 
corresponding component in the crack tip equations must be discarded so that the 
number of equations and the number of unknowns match. 


5.3 Mathematical Formulation and Numerical Method 

The problem of a half-plane weakened by a system of equidistant parallel 
surface cracks with alternative lengths a, and a 2 can still be formulated in terms of a 
system of singular integral equations. The stress equilibrium conditions may be written 
as follows: 


<rA°,y) = A(0 ACO 'Zz( t > b+2nb >y)dt 

** ft n-~<o ** ft n=-oo 


(5-5) 


= ^ I +<*.00 


for 0 < y < a, 


and 


cr. 


(i.y) = r"f D l (t)t,g(‘’ b + 2nb -y'>‘“ + f£ D 1 (t)J j g(,t,2nb,y)dt 

* ft rt=-® ** ft n=-°o 


E (<h 




(5-6) 


for 0 <y <a 2 


where cr, (_y), cr 2 (y) are cohesive stresses at corresponding crack surfaces. The 

function D x ( t ), D 2 ( t ) represent dislocation densities and the kernel g is given by 

(t n _ 2 (y + 1) (y + 3t)[(y + t) 2 -x 2 ] 4 ty(y + t)[(y + 1) 2 - 3x 2 ] 

gK ’ ,y) (y + t) 2 + x 2 [(y + 1 ) 2 + x 2 ] 2 [(y + t) 2 + x 2 ] 3 

(y-t)[(y-t) 2 +3x 2 ] 

[ O '-') 2 + x 2 ] 2 


(5-7) 
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Since g is an even function of x, the series in (5-5) and (5-6) may be written as 

Z g(t, 2nb,y) = g(t, 0 ,y) + g(t, 2 nb,y) (5-8) 


n=l 


Y J g(t,b + 2nb,y) = 2 < £ j g(t,(2n-\)b,y) 


(5-9) 


n= 1 


The above series can be summed in closed form by using same technique in chapter HI 
and written as 


2±g(,,2r.b,y) = zpL + £ coth ^±£- 0 > *2Q± csch > „Z±1 

ti ( y + t ) 3 b 2b 4 b 2 2b 

tyi? y + t . y + t 

+ , csch it- coth7r- 

2b 3 2b 2b 


(5-10) 
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1 - it — — -coth it — — - + n r 2 
b 2b 


2b ) 2b 


2i ^ (^(2 ^ -l)^)=|.anh^ + ^±|^ S ech» ^ 

0 Z0 40 ZO 


n=l 


(y^ 3 
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;r 

46 


se ch ! ^ 
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tanh7r 


y + r 


. , y-f y-t 

4 tanh ;r- + it- 

2b b 


2b 

2 


(5-11) 


sech 2 n — — - 
2b 


It is further convenient to represent the kernel function as 
g{t,Q,y) + 2'Y j g(t,2nb,y ) = ~g, 

n=l O 


itt /rjA 

T’T J 


(5-12) 

(5-13) 


where the function g, , g 2 can be written as 
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g. (t>y) = — - — + 2 coth(^ + t)-(y + 3t) csch 2 (y + t) + 4ty csch 2 (y + /) coth(>> + /) 

t-y 

+ — - — [l -2 (y- 1) coth(y - /) + (y - t) 2 csch 2 ( y - f )] 
y-t L 1 


( 5 - 14 ) 


g 2 (t,y) = 4t + 2 tanh(j + 1) - (y + 3/) tanh 2 (y + t)- 4ty tanh(y + f ) (1 - tanh 2 ( y + /)} 
-2 tanh(_y - /) + (y - 1) tanh 2 (y - 1) 


( 5 - 15 ) 


The crack opening displacements can be represented by the dislocation density D(y) 
[Keer et al. 1975] as 


w i(y) = £ D \(0 dt , 0 <_y <a, 

w 2 O') = f* D z(O dt , 0 <y<a 2 

Thus cohesive stresses at crack faces then can be expressed as 


( 5 - 16 ) 


<rx(y) = fA l — — r D x {t)dt) , y <y<a, 
w c J y ' 

°2 (y) = ft {1 — - r A (')<*} > y <>y<.a 2 

w J y 1 


( 5 - 17 ) 


Because function D x ,D 2 are singular at crack tips, it is now convenient to introduce 
new unknown function C, , C 2 : 


i 



i 



( 5 - 18 ) 


which are smooth functions. Equations (5-5), (5-6) are thus transformed into the 
following equivalent integral equations: 
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Next, extending the interval into (-1, 1) and employing the same technique used in 
chapter III we can convert the singular integral equations into discrete forms as 


n^2n + \ 2b V 2b 2b J 4;r“7 2w + l 2b 


E £ 

4 


= - //, / [^3 + /, { 1 - 1 v , B, (S, )} , 


2b ’ 2b J 

k = \,2,....,n l 

(5-24) 


E V * c ( gM na 2 X l \ I E V £ P _ fejy 1 E£i 

4^~^2n + l 1 ‘ 2b ^ 2 (. 2b ’ 26 J 4n% 2n + l ^ j) 2b g \ 2b ' 2b 


= ~YfJ \^ L ) +/, (1 - v, S 2 (* )}, 

ydj w c ti 


l=\,2, -,rh 


(5-25) 

where v ( are integral coefficients corresponding integral point i (because of the integral 
interval is not (-1, 1) we can't use Erdogan's integration formula for this part). Where 


s, = cosl 


2 / — 1 
v 4/r,+2 


7T , X, = COd 


r _kx_} 
v 2n i + 1 y 




s = cosf — — — n 
1 \^ 2 +2 ) 


( In 
x, = cosj^ 


2/7 2 +l y 


(5-26) 


let 


a \ a 2 A d 

e, =— , e, = — , A = — , 

' b 2 b b 

ne. 


t\ t = ne\Si, t 2j = ne 2 s i , y ]k =— L x k , y 2l = 


_ ne i 


(5-27) 


and introduce 
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4ta)=- 


Msj)=- 


E 1 

4 ( 212 , 4 - 1 )/ 
E 1 


5 ,( 5 ,), / = 




(5-28) 


4(2w 2 + 1) / 

then equation (5-24) and (5-25) may be written as 

Z A (*f )# . ('u .>>1* ; *1 ) + X ^ ( S J ; e 2 ) = r/ — 1 - 2a, * (2/1, 4-l)£v ( 4(s,), 

/=1 i=l v ^ y 7=1 

^ = 1,2, n, 

(5-29) 

X 4 ( 5 . V*2 ('ll > ; e i ) + X ^ ( s y ;e 2 ) = yf -^-7 2i j-l-2a 2 *(2/i 2 +1)5 ]v i . 4(5 i .), 

/ = 1,2,.. ..,nj 

(5-30) 


i=i 


i=i 


where 


H x (t,y;e) = — {4/ 4- 2 coth(>> 4- 1) - (y 4- 3t) coth 2 ( y + 1 ) 

A* 

44 ty coth(_v 4- /)[coth 2 (y 4- /) - 1] - 2 coth^ -t) + (y-t) coth 2 ( y - /)} 

(5-31) 

and 

jrp 

H 2 ( t , y\ e) = — (4/ 4- 2 tanh(_y 4-/)-(_y4-3/) tanh 2 (y 4- 1) 

4-4 ty tanh(j> 4- 1)[ tanh 2 (y 4- /) - 1] - 2 tanh(y - 0 4- {y - 1) tanh 2 (y~t)}. 

(5-32) 

and in deriving the last terms of these two equations the expressions for w c and /„ have 
been taken into account of. 


Also it is convenient to write the stress equilibrium condition in the matrix form as: 

[H]{A}+a 1 *[V](A} = {F} (5-33) 


where 
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”tt tt 


V, 0 ■ 

r 

[H] = 

“ll 

h 21 

. [V] = 

0 fLV, 

L J 

. M=r| 


Note that only at a, = a 2 (or e, =e 2 ) and n x = 



H u = H 12 = H 2I , V, = V 2 , f , = f 2 (5-35) 

In the cohesive crack model the stress intensity factor K is always zero. This crack tip 
condition is equivalent to 

4(1) = 0 (5-36) 

and 4(1) = 0 (5-37) 

This is our second condition that should be satisfied. 4 ( 1 ) a °d 4 > 0 ) can be 
obtained by extrapolating the A(sj). The equations (5-33) and one of (5-36), (5-37) can 
be arranged as a system of equations from which one can solve for A(sj) and the 
loading parameter y . And rest one of equation of (5-36), (5-37) can be used as a 
checking condition. 


5.4 Computational Procedure and Observation 

Through these two conditions or four governing equations we can find some 
very interesting phenomena and fracture behavior. 

To solve the problem, the following procedure are used: 

(1) give value of space b* and load depth d*; 

(2) Select a value of crack one length a,* ; 
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(3) For a guessed value of crack two length a 2 * solve the Eqs.(5-33) and (5- 
36) for {A} and the loading parameter y. 

(4) Adjusting a 2 * to satisfy the equation (5-37) based on the solution obtained 
in step (3), go back to (3). Repeat this loop until all the equations are satisfied. 

(5) Choose another value of a* and go back step (3). 

We did check the crack opening displacement of crack two at every step in the 
unloading condition and see if it gradually close when crack one still keep growth. 

In practice we use the same number of integral points on each crack surfaces, 
thus n x -n 2 =m. Also the size of the system of equations (5-33) is denoted as n, 
which is equal to 2m. The initial results are plotted in Fig. 5-3 The following phenomena 
are observed: 

(1) . In the range from the origin to point Cbif (a) cracks grow stably at the same 
length with increasing of load and (b) no crack closure is observed. The solution (a) is 
checked by setting cracks at different lengths. They come out to be the same and 
satisfy all four conditions, (b) is checked by setting crack two unloaded. But the crack 
opening displacement at crack two still keep growth after the unloading. It means that 
the crack two cannot stop growing at this range. 

(2) . After point B the crack growth has at least the two possible paths: top one 
from point C*,/ to C max in which cracks grow in same lengths and bottom one from Cbif 
to C m i in which one of cracks, say, crack two, stop growing at point Cbif and start to 
gradually close up thereafter, while the crack one continues to grow (possibly at a 
faster rate). It is a possible that crack two unloading at any point between C*,y and 
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Cmax The upper curve corresponds to the condition that the applied load is maximized 
when the crack length is equal to the given value for a given crack spacing. Obviously 
the load required for the unloading path right at bifurcation point (bottom one) is less 
than anyone else. One can easily understand why one crack unloading is more likely to 
happen in the real situation. 

(3). When we solve the system of equation (5-33) we actually first triangulate 
m as ram, where [Q] is lower triangle matrix and [P] is upper triangle matrix 
with Pjj =1. So one can obtain the determinant of [H] by multiplying all Q u . Instead 

of checking det.[H] one only need to check Qjj. It is observed that all Q u are positive 

before point Cb,f, Qnn is almost zero at point C bi f then changes to negative sign after 
Point Cbif. That means det.[H] becoming negative. From C bi f to Cmax the negative Qjj 
moves up from bottom position. When reaching point Cmax there are two Qjj becoming 
negative but det[H] change back to the positive sign. 

On the contrast, det[H]b (take unloading equation (5-4) into account) is 
positive between C*,y and C„i until reaching C„i where Qmm starts become negative. 

Mathematically there must be a some equation, which can decide the position 
Cbif, Cmax and C m! where det.[H] change its sign. This work will greatly simplify the 
fracture analysis. We no longer need step-by-step to search these critical points. 

5.5 Rate Equation of CCM 

To study these solution behavior, one needs to study the rate form of the basic 
equations. For developing the rate equation of CCM in the initiation problem with two 
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cracks the stiffness influence functions will be utilized in our basic conditions. The 
stress equilibrium conditions at crack faces can be written in terms of stiffness function 
as 


- F(x, ) + o< w, ( x, )) = JT S x , ( x, , V )w, ( x, ' )<&, ' + P S l2 (x, , x, ' )w 2 (x, ' )dx 2 ' ' 

Jo a Jo ( 5 - 38 ) 

- F(x 2 ) + o{w 2 (x, )) = £ 5 j, (Xj , X, ■ ' )w, (x, ' )<&, ' + JJ 1 ^ 22 (X2 , X3 ' )w 2 (Xj ’ ’ 

where x, (i = 1,2) are the coordinates measured along crack one and crack two. The 
stiffness functions S\y(x it x } ) ( force at x, generated by a unit displacement at x ; .) is 
symmetric with respect to i and j due to the assumption of linear elasticity in the bulk 
material. The above conditions can be simplified by using the summation convention as 

G- = F(x f ) - a{w i (x,)) + £ J S IJ (x,,x)w j (x)cix = 0 (5-39) 


where i =1, 2. The summation convention in this section only apply to the repeated 
subscript j. 

The total stress intensity factor in each crack tip must be zero, which can be 
expressed as 


K x = £ k u (x)w, (x)dx + £' k n (x)w 2 ( x)dx = 0 
k 2i OK (*)<& + k 2i ( x ) w 2 ( x )ckc = 0 


( 5 - 40 ) 


where & l7 (x y ) = stress intensity factor at the tip of crack i due to a unit displacement at 
the position Xj. Also crack tip conditions can be simplified in same way as 

^ = r ' k i j {x)w j {x)dx = Q 


( 5 - 41 ) 
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Equation (5-39) and (5-41) are the basic equations of CCM with two cracks. If 
there are more than two cracks, one can simply extend the range of index j. If the range 
of j is from 1 to n, then /' = l,2,...,2n, so there are just enough number of equations to 
solve the unknown crack opening displacement Wj and crack length aj , It is useful to 
call q t = (Wi , a, ) the state variables, for the system is totally determined once the state 
variable is specified. 

The rate equations can be obtained by considering the basic unknown variables 
to be the function of time t, which is only for the purpose of keeping the sequences of 
the system development. When a, changes the stiffness functions also change. 

The derivatives of the stiffness functions with crack length a, can be expressed 
in terms of stress intensity factors as 

( 5 -«) 

M 

No summation over m is implied. This relation is given, for instance, by Okamura 
(1975), and (5-42) is simply a generalization to the case of multiple cracks. 

The time derivative of (5-39) can be expressed as 

+ J7 ■S,( x i’ x )”V( x ) £& 

. Z; (x . x) ( 5 - 43 ) 

+ ~ \ ~ M 'j( x ) dx + a j S,J ( x , » a j ( a j ) = 0 

where a dot denotes derivative with respect to time t. & denotes the derivative of 
softening function a with respect to its argument w. However, at crack tip 
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displacement Wj is zero, thus the last term can be dropped. The next last term can be 
changed due to (5-42) and (5-41) as 


<[ J ^ SlJ £ , ’ X ^ w J (x)dx = k mj (x)w j (x)dx = —djc^x,)^ ( 5 - 44 ) 


This term is also equal to zero due to K m = 0 
Thus, equation (5-43) reduces to the form 

dG ; 


-^ = -a , w i (x i )~ S tJ (x . , x)w } . (x)dx = 0 


( 5 - 45 ) 


which is still independent of a m . If there are unloading cracks, the rate equations should 
be modified accordingly. In other words, & should be calculated according to 
unloading curve. 

The time derivative of (5-41) can be expressed as 

= £ k ,j (*)”', (x)dx + Kj ( x >j (x)<& = 0 ( 5 - 46 ) 


In calculating the last term, we deliberately avoid expanding the derivative into the sum 
of the derivative with respect to the integral limit and the derivative of k,j with respect 
to a„, because each term will be unbounded. However, since we expect the final result 
to be bounded these singular terms must cancel each other to yield a bounded value. 

Equations (5-45) and (5-46) are the rate equations of CCM for the thermal 
loading problem. They depend on the state variable , but is linear in the rate of state 
variables. When the system is normal, there is only zero solution to the rate equation. 
When a system is critical, there can be a non-zero solution to the rate equation. In this 
way, one can study the behavior of CCM by study the behavior of its rate equation. 
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One can find the rate equations (5-45) in the matrix form as 

[H]{A) +a*[V]{A} = {0} (5-47) 

To solve it for a non-zero solution is a difficult task. However, the problem is 
simplified if (5-47) is viewed as an eigenvalue problem. In particular,, we are interested 
only in the smallest eigenvalue of a[. 


5.6 Maximum Load and Bifurcation 

A symmetric solution, in which the two cracks are of equal length and same 
crack displacement, was plotted in Fig.5-3 as the curve on top. The system is at critical 
condition when the applied load is at its maximum, and the corresponding crack length 
is denoted as Cmax, which is smallest eigenvalue that corresponds to the eigenvector 
{A,} = {A 2 } (equivalent to w, = w 2 ). For different crack spacing these crack lengths, at 
which the applied load is at its maximum, compose a right curve in Fig 5-4. 

The maximum load can be calculated directly using the eigenvalue solution. It is 
obtained by multiplying the stress equilibrium condition (5-39) with eigenfunction and 
integral it along the crack surface as 


Of 

I rff ( X /d) ~ ft ( 1 - — ) + £ S tJ (x, , x)w J . (x)cix 
0 L W c 


w,. (*,)<&,. = 0 


(5-48) 


By changing the integral order and using the rate equation (5-46) to cancel relevant 
two terms, we find that the loading parameter can be computed from the expression 


r = 


\°™(x)dx 

^M'{x)f{xld)dx 


(5-49) 
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In deriving this equation, we have utilized the critical condition (5-46) and the 
symmetric property of the stiffness functions, as well as the symmetry condition of the 
eigenfunction. This solution is useful if one need to know the maximum load as a 
function of crack length. 

The left curve in Fig. 5-4 corresponds to the bifurcation point for a given crack 
spacing, and the corresponding crack length is denoted as C*,/ in Fig.5-3. This curve 
corresponds to the eigenvector {A,} = -{A 2 } (equivalent to w, = -vi» 2 ), so the 
corresponding load parameter cannot be calculated using a relation similar to (5-49). 
The symmetric solution is the unique solution if or/ is less than the bifurcation length. 
When passing the bifurcation length, one of the cracks may unloading, leaving only the 
other one to grow in response to the loading. Consequently, the solution lost its 
symmetry. The load required to propagate only one crack is usually less than that of a 
symmetric solution. Although both solution are stable , the asymmetric solution is more 
likely to be actually followed. 

5.7 Maximum Load with One Crack Unloading 

The unloading stress-displacement equation (5-4) need to be used in crack 
unloading condition. The solution procedure must be modified because zero K 
condition for unloading crack must be discarded. The load parameter as a function of 
the active crack length is shown in Fig.5-3 as the curve in the bottom. 

With one of the cracks unloading and the other growing, load reaches a 
maximum value, which corresponds to yield another critical condition. The critical 
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condition can again be characterized by the rate equation (5-45) or its matrix form (5- 
47), in an expanded form, as 


[H 11 ]{A,} + [H„]{A 1 } = - a ;[v 1 ]{A,} 

(5-50) 

[H 11 ]{A,} + [H n ]{A 1 } = a;[V ! ]{A 1 } 

(5-51) 


where [v%] is modified according the unloading stress-displacement equation (5- 
4).This system equation cannot be simplified into a proper linear eigenvalue problem, 
because there is a negative sign in right hand side. But one can eliminate {A z } from (5- 

51) 

{A 1 } = -([H K ]- a ,-[V ! ])" l [H !l ]{A 1 } = -[Q]{A,} (5-52) 

[Q] can be solved from the following equation 

([H B ]-af[V J ])[Q]=[H„]. (5-53) 

substituting the (5-52) into (5-50), a well-posed eigenvalue equation for solving {A,} 
can be obtained as 

([H, ! ][Q]-[H„]){A,} = < I ;[V 1 ]{A i } (5-54) 

The solution strategy is: (1) assume initial a*; (2) calculate [Q] from (5-53); (3) solve 
the eigenvalue problem through (5-54), and corrected value of a* obtained; (4) use this 
new a* and repeat step (2) and step (3) until the value a\ is converged. This procedure 
usually takes only very few steps, typically 2 or 3 for 4 digit accuracy in a*. Since the 
singularity condition is mainly determined by crack one and the influence of crack two 
on crack one is not very substantial. 
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The smallest eigenvalue of the singularity condition yields yet another critical 


crack length for crack one, which corresponds C m j in Fig. 5-3 and the curve plotted in 
the middle of Fig. 5-4. It is interesting to note that this curve is virtually the same as the 
right curve if its crack spacing b* is multiplied by a factor of 2. In fact, the crack 
spacing should be doubled when the crack one keeps growing and the crack two 
unloads at early stage. I would explain it later. The corresponding maximum load can 
also be obtained by multiplying the eigenfunction to the equilibrium equations and 
integrating along the crack surface 


['M*)dx 

J 0 °' w. (*)/(* / d)dx + £* w 2 (x)f(x / d)dx 


(5-55) 


In above section, we have assumed that crack two starts to unloading right at 
the point of bifurcation. Although this is most likely the case, it does not have to be this 
way. Actually, a 2 can be any value between C^f and C max . In this sense, the solution is 
not unique once the crack lengths exceed the bifurcation point. This is very similar to 
the case of the buckling of an elastoplastic column (Bazant and Cedolin, 1991): once 
the applied load exceeds the tangent modulus load, there are more than one solutions 
bifurcating from the original symmetric solution. 


5.8 Post-critical Behavior and Lower Crack Spacing Limit 

We are also interested in the crack growth after these critical points. From 
Fig.5-5a the view is very clear that, there will be a jump in crack growth if every 
another crack is unloaded at the bifurcation point. Note that the crack opens fully at 
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a, * = 0.734 , which is half way during the jumping, therefor it cannot stop at that point. 
The crack length at the post-initiation state should be measured from the surface to the 
point Cm, where the load once again becomes the same value as the peak load. Before 
this point crack would not stop and after this point the crack would -not grow without 
further increase of load level. The crack mouth opening displacement of crack two is 
reduced to below 0.1 while that of crack one is increased to 2, see Fig.5-5(b). When 
crack two is compared with crack one, both crack length and crack opening 
displacement are very small. This means that practically the cracking spacing is 
doubled. The equal length cracks (upper curve) will grow slowly and require much 
higher level of load to initiate the cracks will not actually happen. 

When the crack spacing is small, like the previous example, the cracks become 
highly interactive, and some cracks stop while the others grow at a much faster rate. 
When the cracks are far apart, say b*>20, their interaction is weak and , therefore, they 
grow at equal rates. A question therefore naturally arises: does the crack spacing have 
a lower limit? The lower limit is defined as the least crack spacing the cracks can 
initiate. In other words, even denser cracks may exist in initial stage, some of them will 
stop in the process of crack initiation and thus should not be counted when measuring 
the spacing in the final stage of crack initiation. 

The lower limit of crack spacing is found to be 2, as can be seen from Fig. 5-4. 
This lower limit for quasi-brittle material in cohesive crack model is the same as the one 
for brittle material discussed in the previous chapter. The reason is as follows. The 
process zone length, although varies with structure size, has an upper limit. When a 
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crack length is very long and the relative process zone becomes vanishing small. As a 
result, the behavior of the cohesive crack model converges to that described by linear 
elastic fracture mechanics when the crack length becomes very long. For the lower 
limit, the crack depth approaches infinity, that is why both theories haye the same lower 
limit value for crack spacing. 

5.9 Concluding Remarks 

The crack initiation is a complicated problem. Although we have established a 
theory to determine initial crack spacing based on linear fracture mechanics, such a 
theory cannot be applied when material is not perfectly brittle. In real world, there is 
no perfectly brittle material. Cohesive crack model is the simplest possible nonlinear 
fracture mechanics available to study this problem. With cohesive crack model, the 
crack tip is replaced by a nonlinear process zone and the concept of fracture is 
extended. With such a model it is possible to treat the initiation problem at a more 
realistic level. 

It is observed that there is always a considerable amount of randomness 
associated with initial crack spacing. In reinforced concrete, the crack spacing is 
described, not by a determined value, but by a lower limit. In the light of cohesive 
crack, the initial development of the nonlinear process zone obviously will start from 
the weakest point of material. Those weakest point of the material are, usually called 
the defects, are randomly distributed. Such a point of view is adopted by Weibull 
(1951) in his analysis of statistical behavior of material strength. However, as we have 
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known from our preliminary study, not all defects will develop into macrocracks. The 
geometry as well as the loading process has an effect on which defect will become 
activated. Our main goal is to understand the basic mechanisms of this selection 
process and find a way to model these mechanisms. 

The crack initiation described by the cohesive crack model with alternative 
crack length is much more complicated than the one for brittle materials. It is a gradual 
process that may involve propagation as well as crack closure. In the example 
problem, both bifurcation and stability limit can happen before the crack mouth opening 
exceeds the threshold value Based on the rate equations of the cohesive crack 
model, mathematically it is shown that the critical • conditions of the cohesive crack 
model can be determined by considering only the homogeneous rate equation of the 
stress equilibrium equation. 

The external loading providing the energy for crack initiation includes two 
parameters. They adjust each other: the longer the penetration depth d, the less the 
maximum loading stress y. Because the homogeneous rate equation does not include 
the external loading, crack lengths of bifurcation, critical points determined by the rate 
equation depend only on the crack spacing, except when there is an unloading crack 
(but the effect due to the external loading in the unloading crack has negligible 
influence on the crack spacing, as has been mentioned before). 
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Fig. 5-1 (a) Cracking of quasi-brittle materials 


a 
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Fig. 5-1 (b) The cohesive crack model 



Fig. 5-2 (a) Geometry difinition of parallel crack system 
with alternative lengths; 


(b) Unit cell of width 2b. 
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Crack Length a^ 


Fig. 5-3 Possible crack initiation paths 



Crack spacing b* 



Fig. 5-4 Bifurcation and critical curves 



Crack Mouth Opening w/w Loading Parameter 


b* = 2.5, d* = 2.5 
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Fig. 5-5 Post— critical behavior ot crack initiation 




CHAPTER VI 


CONCLUSIONS AND FUTURE RESEARCH 

The study of crack initiation from a smooth surface is still at its infant stage, 
many problems are still not fully understood and require further researches. However, 
now it has been established in this study that in order to describe the initiation, the 
strength criterion and the fracture mechanics theory must be combined. Strength 
theory can only mark the starting point of cracking, while the fracture mechanics 
determines what would happen after the cracking is started. In the general situation, 
the material under consideration always have defects to begin with, so the strength of 
the material varies from point to point, and cracks of small sizes will start where the 
strength criterion is first reached. In this sense, the crack initiation is a random 
phenomena. However, not every small crack develops into a crack of macroscopic 
sizes, most cracks will be unloaded and yield themselves to a few cracks of larger sizes. 
In other words, , the structure appears to have the ability to select, from a large number 
of cracks of microscopic sizes, a few cracks that will eventually grow into macroscopic 
sizes. Although the overall initiation process is random, this structural effect is 
deterministic in nature and is the main reason why its average behavior can be studied 
by a deterministic approach. 
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If the material can be modeled as linear elastic and perfectly brittle, the crack 
initiation can be summarized as three conditions between the two states of the initiation 
process. The first condition is the strength criterion, which determines the stress level 
at which the pre-initiation state. This stress level is assumed to be kept constant during 
the whole the initiation process. The second is the condition that the Griffith condition 
be satisfied at the post-initiation state. This condition is important not only because it 
marks the end of the initiation process, but also because it determines the total amount 
of the elastic strain energy the structure can release. The third condition is required 
due to the fact that the crack initiation is not an equilibrium state. Cracks are often 
occur suddenly, indicating the process is more like a jump than gradual growth. In 
fact, because the energy release rate is a linear function of crack length if the length is 
very small, the Griffth condition cannot be satisfied during the early stage of the crack 
initiation. That is why it is essential to supplement the condition that the total energy 
released during the process of crack initiation must be converted into surface energy. 

The study of crack initiation from the surface of an elastic half plane is related 
to the crack initiation from the surface of a large structure, such as the top surface of a 
large concrete block. The solution can be related to the spacing and the depth of 
cracks that are some distance away from the side wall. It is interesting that the crack 
spacing has a lower limit of 2 1 0 . The same lower limit of crack spacing can be obtained 
by using just the first and the third condition. This result should be explored further. 
The case of infinite crack spacing warrants further attention. According to our 
assumption, the crack will not propagate if the loading depth is very short, even though 
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the stress level may exceed the tensile strength. Further experimental study is needed 
to clarify the situation. 

The problem of crack initiation in the pavement is interesting not only from 
mechanics point of view, but also from the practical point of view. Cracking is one of 
the main factors that cause damage to the highway pavement. The study represents the 
first step to a better understanding of the behavior of cracking in the pavement. For 
Portland concrete pavement, joints must be provided in the pavement so that random 
cracking can be controlled. The theoretical result may help engineers to design the 
joint placement on a physically sound basis. 

From the mechanics point of view, pavement is a non-homogeneous structure: 
the material properties of the top layer is quite different from its subbase. Therefore, 
the half-space solution cannot be directly applied. However, since the typical crack 
spacing is large compared to the pavement thickness, the plate assumption can be 
advantageously be used to simplify the analysis. The nonlinear process zone is 
important because the zone size is not small compared to the pavement thickness. As 
an approximation, it is assumed that in the third condition, the fracture energy should 
be determined by the work-of-fracture method, while in the second condition, the 
effective fracture energy. The modification is necessary because the effective fracture 
energy is usually significantly different from the fracture energy calibrated by the work- 
of-fracture method for pavement material. 

This brings us to the discussion of crack initiation by the cohesive crack model. 
Concrete and asphalt concrete are classified as quasibrittle materials, meaning that 
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cracking is the main failure mechanism. In the cohesive crack model, the crack is a 
gradual process. Once the tensile stress reaches the tensile strength of the material, the 
process zone, where bridging stress decreases with increasing crack open displacement, 
starts to elongate until some limit size is reaches which corresponds to the maximum 
loading value. The first condition remains the same if one assumes that the pre- 
initiation state can be marked by the starting of the process zone. The second 
condition is essential unchanged, except now it is in the form of the condition of zero 
total stress intensity factor at the tip of process zone. The third condition can also be 
translated into the language of the cohesive crack model if we could define the total 
energy of the cohesive crack model. However, the total energy cannot be meaningfully 
introduced if there is crack closing. This is perhaps the main reason why more study is 
needed along this direction. Nevertheless, the cohesive crack model can successfully 
be used to identify the moment when uniform cracks growing bifurcate to only every 
another crack grow while the rest stop and unloading. Such phenomenon is believed to 
be important to fully understand the crack initiation process with more realistic detail. 
When the loading depth is increased, it is seen that both cohesive crack model and the 
theory of linear elastic fracture mechanics yields the same lower limit for crack spacing, 
no matter whether unloading is considered or not. 

For further research directions, we may note that the crack initiation theory has 
not been established, as yet, for the cohesive crack model. The difficult lies in the fact 
that the very concept of initiation becomes blurred once the cohesive crack model is 
introduced. Because the material is still stable when the process zone is just started, 
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the stress level must be raised further to make the cracks to grow. In other words, 
there is no jump of crack propagation before the maximum load is reached. Since there 
is no stage during which the incremental energy balance condition is violated, the third 
condition of the crack initiation theory may not be an independent condition, at least 
from the physical point of view. 

In addition to the problems that related to the cohesive crack model, the theory 
still needs to be further generalized to cover more complex situations, among them the 
following can be mentioned: 

1) The crack initiation theory of composite material. One of the main energy 
absorption mechanism is multiple lateral cracking along the fiber-matrix interface. 
Dense multiple cracking is a desirable phenomena in certain applications. It seems that 
the research idea developed here can also find application in that area. The main 
modification would be the consideration of the shear force that is transferred from fiber 
to matrix or vice versa. The crack initiation problem in reinforced concrete belongs to 
the same category. 

2) In this study we have only considered the solution of one dimensional crack 
initiation pattern. The pattern in two-dimensional plane is known to hexagonal, but the 
proper equation has yet been set up for this interesting problem. 

3) Another generalization is the problem of non-uniform solution. We have 
been able to obtain the uniform solution as the average behavior of cracking pattern, 
but it will be a big step further if we can also know something about the non-uniformity 
of the crack pattern, and understand what we can do to control the situation. 
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4) The problems discussed in this study is defined on a surface extending 
infinitely in both direction. If there is boundary that ends the extension of the free 
surface, then the crack initiation theory developed here must be modified, because the 
crack pattern may not be uniform any more. This problem is thus also connected to the 
previous problem. For a non-uniform pattern of cracks, the parameters that 
characterize the pattern is more than two, so the theory developed herein is not 
sufficient to determine any crack pattern that needs more than two parameters to 
define. This seems to be the most urgent problem that must be solved. 

5) At present research stage only mode-I crack is considered. For mixed-mode 
cracking (for example, half-plane with a system of parallel incline cracks) the shear 
cohesion should also be included. The stress-displacement relation for the shear 
cohesion may have a similar form as the one for the tensile cohesive stress. The number 
of the basic equations are doubled but the solution procedure remains the same. 
However, since shear deformation is often associated with friction, which consumes the 
released energy, the energy conservation equation must be modified accordingly. 
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APPENDIX 1 


THE STRESS INTENSITY FACTOR 


We note that, except for the first term in gj (given by equation (3-14)) function 
g! is bounded. (The last term in g l is bounded that can be proved by expansion of 
hyperbolic functions in power series.) Also the unknown function C is a bounded 
function. Once the unknown function is solved, the stress intensity factor can be 
calculated as 


y->a 




: lim — , 
4/r 

COO 



( 2 \ 

J>) 


{a 2 -t 2 j 


1/2 


dt 

t-y 


=dt + r 

! JO 


• C(t)-C(y) 


(t-yyja'-t* 


dt 


(1-1) 


If we prove any term in the brace is bounded the corresponding product in K 
will be equal to zero because lirn yj(y - a) = 0 . 


y—>a 


Now we prove the last term in brace is bounded function. 

C(t) is a smooth function that will satisfy the Holder condition (Muskhelishvili, 
1992, page. 11) : 

\C{t)-C(y)\<M\t-y[ (1-2) 

where M and p. are positive constants. 
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By using this condition the absolute value of the last term in the brace can be 
found as 


I 


CU)-C(y) 


dt<M r r f - , — (1-3) 

J o \t-y\ J a 2 -t 2 


\(t-y)yla 2 -t 2 

Because the exponents (l-/i)<l, this generalized integral is finite (bounded). As a 
result, the corresponding term in K is equal to zero and only the first term in the brace 
need further consideration. 


The first term in the brace can be expressed as 
C(y ) . . r° dt 


The last term in equation (1-4) obviously is a bounded function. For the first term in 
equation (1-4), we first normalize the integral interval as (r = t/a): 

C(y) r> dr 


op 

a (r - y / a)Jl-r 2 


(1-5) 


When y -» a + according to Delale and Erdogan (1977), the integral above is equal to 


COO 


— TCI 




+ *00 


(1-6) 


where function R is bounded at a + . So the final result of AT is 


K = lim — ^2 7t(y - a)aC(y) 

y-*a* 


-IK 




= --.JanC(a) 


(1-7) 


APPENDIX 2 


ADDITIONAL COMPLIANCE FUNCTIONS 

The additional compliance functions can be calculated according the formulas given 
by (4-42) to (4-44) using numerical integration method. Where the stress intensity factors 
are defined by (4-25) and (4-26). Specifically, we use Simpson's rule to perform the 
quadrature with the integral domain [0,1] subdivided into 200 intervals. The obtained 
numerical results are used as a basis to construct approximate expressions for the 
compliance functions. 

The nondimensional compliance function for tension, given by the handbook, is 

4W = 7 ^-j[0.99-x(1-x)(1.3-1.2x + 0.7x 2 )] (2-1) 

(1-x) 

which is very close to our numerical integration result. Using the following form of the 
expression, we obtain the optimum fit as 

Aw(x) = — [0.988 -x(l-x)(l. 504- 1.654x + 0. 920 * 2 )] ( 2 - 2 ) 

(1-x) 

Although visually indistinguishable, (2-2) is more accurate than (2-1). 

For bending, the compliance function given by the handbook is 

= — - 7 (5.93-x(19.69-37.14x + 35.84x 2 + 13.12x 3 )) (2-3) 

(1-x) 
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which correlates reasonably well with the result of numerical integration. However, the 
following expression yields a far better approximation to the numerical results: 


tW - — j[0.988+*(l-x)(0.413 + 0.282x)] (2-4) 

(1-x) (l + 2x) 

especially in the range where x is close to 1 . 

For the mixed compliance function, there is no expression listed in the handbook. 
Based on the result of numerical integration of the stress intensity factors, we obtain the 
following expression: 


Aaw (*) 


24; 


(l-x) 2 (l + 2x) 


[0. 988 - x(l - x)(0. 578 - 0. 689x + 0. 196 x 2 )] (2-5) 


In the actual calculation, we use our own expressions (2-2), (2-4) and (2-5). The forms of 
these expressions ensure the correct asymptotic limits for at both x -> 0 + and x -» 1“ . 


